Program VALUEFUNCTION

USE PARAMETERS
USE Functions_libs
USE My_ToolS
USE Interpolate_func
use omp_lib
Implicit NONE

Integer :: jo_lagged, joption_pay, jfix_par, jht_par, jpay, jhc_ret, jmc, jtrpay, jtype4, jtype4_p, jl_a, jm_a, ju_a,  jl_hc, jm_hc, ju_hc, jhc, jexp, jDI, jDI_p, jos, jos_p, jMEcat_p, jMEcat, jExperimentNumber, jspouse, counter, max_hk_age , jzhc_p,jzhc, jmar, jmar_p, jinc, jtype3,jtype3_p, mhk,jzt,jzt_p, kk, jas,jas_p,jsav,je,jh,jh_p, jr,jr_p,ji, jage,j, jht, jyears, jfix, jtype1,jtype2,jtype2_p,js,jdp,jdu,jx,jo,jo_p,  jins, jhrs , jx_age, jsav_upper_ret  , jsav_upper    , jinv_upper  , TR_indic , jhk


Real (Kind=dd), dimension (n_as,n_type1,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime  
Real (Kind=dd), dimension (n_as,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime_simple   
Real (Kind=dd), dimension (n_as,n_type3,n_hc_grid,n_ret-1,n_zt,n_o,n_type4)::  EXP_L_prime_fix  

Real (Kind=dd), dimension (n_o,n_type4,n_zt,n_as,n_type3,n_type1,n_hc_grid,n_ret-1)::  B 

Real (Kind=dd), dimension (n_o,n_MEcat,n_type4,n_zt,n_as,n_type3,n_type2,n_hc_grid,n_ret-1)::  B_1 
Real (Kind=dd), dimension (n_o,n_MEcat,n_type4,n_zt,n_as,n_type3,n_type2,n_hc_grid,n_ret-1)::  B_2 

Real (Kind=dd), dimension (n_trpay,n_o,n_MEcat,n_type4,n_zt,n_as,n_type3,n_type2,n_hc_grid,n_ret-1)::  B_fix 

Real (Kind=dd), dimension (5,n_e,n_ht,n_type3,n_type2,n_MEcat,n_type2  ,n_type3)   :: prob_s_treat_new_64   
Real (Kind=dd), dimension (n_ret:n_age-1, n_e,n_ht,n_type3,n_type2, n_mar, n_MEcat, n_type2,  n_type3, n_mar )   :: prob_s_treat_new_ret

Real (Kind=dd)::  EXP_L_prime_INT, VF, V_death
Real (Kind=dd), dimension (2:3) :: VF_tr
Real (Kind=dd), dimension (n_sav,n_h,n_r,n_o,n_mar,n_zhc):: SAVE_EXP_L_prime_INT

Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_mar):: RET 
Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_ret : n_age,n_mar):: RET_PRIME 
Real (Kind=dd), dimension (n_MEcat,n_h,n_r,n_type2,n_mar):: RET_PRIME_INT

Real (Kind=dd), dimension (n_sav,n_MEcat,n_h,n_r,n_type2,n_mar):: SAVE_RET_PRIME_INT

Real  (Kind=dd) :: Resources, prob_s, prob_d, income_bt, income, base_1, base_2 , Total_tax, actual_earn, earnings_interpol, tax_SS, tax_med, earnings_capacity
Real (Kind=dd) :: last_value, maxB, valret   , check
Real (Kind=dd) :: UTI,x, EXPECTED_RET_PRIME , EXPECTED_B_PRIME_C_FLOOR_NT

Real (Kind=dd)::  slope_hc  , inter_hc    , slope_as, inter_as, HK_temp, CON, TR, ASP , MAX_SAV  , as_opt, TR_opt ,  INT_sav_upper_ret   , INT_sav_upper, EXPECTED_B_PRIME, EXPECTED_B_PRIME_C_FLOOR
Real (Kind=dd), dimension (n_o) ::  EXPECTED_B, jx_p_o
Real (Kind=dd), dimension (n_zhc,n_o) ::   HK_INTERPOL
Real (Kind=dd), dimension (n_e) :: scale_param

Real (Kind=dd), dimension (2) ::  income_pay, Total_tax_pay, base_2_pay, base_1_pay, CON_pay  
Real (Kind=dd), dimension (2) :: MAX_SAV_pay, INT_sav_upper_pay  
Real (Kind=dd), dimension (2:3) :: prob_s_treat  
Real (Kind=dd), dimension (3) :: UTI_trpay, maxB_trpay,  EXPECTED_B_PRIME_trpay, last_value_trpay, V_cf_trpay, V_c_all_trpay
Real (Kind=dd), dimension (n_sav,3) :: SAVE_EXPECTED_B_PRIME_trpay
Integer, dimension (2) :: jsav_upper_pay

Real (Kind=dd), dimension (7) :: x_r
Real (Kind=dd), dimension (14) :: x_h   
REAL (Kind=dd) :: index_h
Real (Kind=dd), dimension (7) :: x_du
Real (Kind=dd), dimension (9) :: x_dp
Real (Kind=dd), dimension (7) :: x_s
Real (Kind=dd), dimension (n_e,n_hc_ret) :: x_ss

Integer, dimension (n_as,n_type1, n_ret, n_o, n_sav) :: AS_LB
Integer, dimension (n_hc_grid, n_e, n_type2, n_h, n_zhc, n_o, n_ret-1) ::  HK_LB  

integer m,n
real(dbl) yt(4), mx1, mx2
 
sav_grid(:)=0.0d0
 do jsav=1, n_sav
 sav_grid(jsav) = dble(jsav)
 end do

do jExperimentNumber = 4, 4
    
 
   

prem_ephi(:,:)=0.0d0 
prem_mcare(:)=0.0d0

open(202, file='Internal_Parameters/parameters_assumed.txt')
read (202,*) interest
read (202,*) taxc
read (202,*) tax_lambda(1) 
read (202,*) tax_tau(1)
read (202,*) tax_lambda(2) 
read (202,*) tax_tau(2)
read (202,*) tau_ss
read (202,*) tau_med
read (202,*) yhat
read (202,*) penalty
read (202,*) prem_ephi(3,1)  
read (202,*) prem_ephi(3,3)  
read (202,*) prem_mcare(1) 
close (202)

if (jExperimentNumber.eq.15) THEN
    taxc=0.065d0
end if 

if (jExperimentNumber.eq.17) THEN
    taxc=0.071d0
end if 


equiv(:,:)=0.0d0
open(278, file='External_Parameters/parameters_fam_size.txt')
do jage=1, n_ret-1 
read (278,*) equiv(1,jage)
equiv(1,jage) = 1.5d0 + (equiv(1,jage) -2.0d0 )* 0.3d0 
end do 

  close (278)

 do jage= n_ret, n_age
   equiv(1,jage) = 1.5d0 
   equiv(2,jage) = 1.5d0 
   equiv(3,jage) = 1.5d0 
 end do
 
 
  
 p_e(:)=0.0d0
open(208, file='External_Parameters/Educ_dist.txt')
do je=1, n_e 
read (208,*) p_e(je)
end do 
close (208)

p_surv(:,:,:, :)=0.0 
open(195, file='External_Parameters/parameters_survival.txt')
 do jh=1, n_h ;  do jmar =1,n_mar ; do jage=26, n_age 
   read (195,*) p_surv(jh,1,jmar,  jage)
 end do   ; end do    ; end do
 close (195)
 
 
   p_surv(:,2,:,  :)=p_surv(:,1,:,  :)
   p_surv(:,3,:,  :)=p_surv(:,1,:,  :)

 
   do jage=2, 25 
   p_surv(:,:,:,  jage) = p_surv(:,:,:,  jage-1) + p_surv(:,:,:,  26)/ 25.0d0
   end do 
   p_surv(:,:,:, :)=1.0d0-p_surv(:,:,:, :) 

    In_mar(:, :)=0.0d0
open(193, file='External_Parameters/Marriage0.txt')
 do jh=1,n_h 
   read (193,*) In_mar(1, jh)
 end do
close (193)

    Tp_mar(:,:,:,:,:,:,:)=0.0d0

open(199, file='External_Parameters/Marriage1.txt')
do jage=1,n_ret-1;   do jo=1,n_o
 read (199,*)  Tp_mar(2,1,1,jage,1,1,jo)
 Tp_mar(2,1,1,jage,1,:,jo)=Tp_mar(2,1,1,jage,1,1,jo)
end do; end do
close (199)

Tp_mar(2,1,2,:,1,:,:) =Tp_mar(2,1,1,:,1,:,:)
Tp_mar(2,1,n_e,:,1,:,:) =Tp_mar(2,1,1,:,1,:,:)

Tp_mar(2,1,:,:,2,:,:) =Tp_mar(2,1,:,:,1,:,:)
Tp_mar(2,1,:,:,n_h,:,:) =Tp_mar(2,1,:,:,1,:,:)

open(200, file='External_Parameters/Marriage2.txt')
 do jage=1,n_ret-1 
 read (200,*)  Tp_mar(1,2,1,jage,1,1,1)  
 Tp_mar(1,2,1,jage,:,:,:)=Tp_mar(1,2,1,jage,1,1,1)
 end do
close (200)
 
Tp_mar(1,1,:,n_ret:n_age,:,:,:) =1.0d0
Tp_mar(2,1,:,n_ret:n_age,:,:,:) =0.0d0

open(155, file='External_Parameters/Marriage3.txt')
do jage=n_ret, n_age-1
 read (155,*)  Tp_mar(1,2,1,jage,1,1,1)
 Tp_mar(1,2,1,jage,:,:,:)=Tp_mar(1,2,1,jage,1,1,1) 
end do
close (155)




Tp_mar(1,1,:,:,:,:,:)=1.0d0-Tp_mar(2,1,:,:,:,:,:)
Tp_mar(2,2,:,:,:,:,:)=1.0d0-Tp_mar(1,2,:,:,:,:,:) 

prob_os(:,:,:,:,:)=0.0d0

open(158, file='External_Parameters/Prob_spouse_work.txt')
 do jage=1,n_ret-1
         do jh=1,n_h
 read (158,*) prob_os(jage,1,jh,2,1) 
 prob_os(jage,1,jh,2,2)=prob_os(jage,1,jh,2,1) 
 prob_os(jage,1,jh,2,3)=prob_os(jage,1,jh,2,1)  
         end do         
end do
close (158)

prob_os(:,:,:,1,:) = 1.0d0 - prob_os(:,:,:,2,:) 

inc_spouse(:,:,:,:,:)=0.0d0 

open(156, file='External_Parameters/Income_spouse1.txt')
 do jage=1,n_ret-1
         do jh=1,n_h
             do jtype4=2,n_type4
 read (156,*)  inc_spouse(jage,1,jh,jtype4,1) 
inc_spouse(jage,1,jh,jtype4,2)=inc_spouse(jage,1,jh,jtype4,1)
inc_spouse(jage,1,jh,jtype4,3)=inc_spouse(jage,1,jh,jtype4,1)
             end do
         end do         
 end do
close (156)


open(157, file='External_Parameters/Income_spouse_ret.txt')    
 read (157,*)  inc_spouse(n_ret,1,1,2,1)     
 inc_spouse(n_ret+1:n_age,1,1,2,1)= inc_spouse(n_ret,1,1,2,1) 
 inc_spouse(n_ret:n_age,1,1,3,1)=inc_spouse(n_ret:n_age,1,1,2,1) 
 inc_spouse(n_ret:n_age,1,2,:,1)=inc_spouse(n_ret:n_age,1,1,:,1) 
 inc_spouse(n_ret:n_age,1,3,:,1)=inc_spouse(n_ret:n_age,1,1,:,1)   
close (157)

inc_spouse(n_ret:n_age,:,:,:,2) = inc_spouse(n_ret:n_age,:,:,:,1)
inc_spouse(n_ret:n_age,:,:,:,3) = inc_spouse(n_ret:n_age,:,:,:,1)

! SPOUSE MEDICAL OOP
oop_spouse(:,:,:,:)=0.0d0  
open(159, file='External_Parameters/OOP_spouse.txt')
         do jage=1, n_ret-1
read (159,*)  oop_spouse(jage,1,3,1) ! has insurance -- own or husband or both
read (159,*)  oop_spouse(jage,1,2,1) ! no insurance
oop_spouse(jage,1,2,3)=oop_spouse(jage,1,3,1) ! insured
oop_spouse(jage,1,2,5)=oop_spouse(jage,1,3,1) ! insured
oop_spouse(jage,1,3,2:5)=oop_spouse(jage,1,3,1) ! insured

oop_spouse(jage,1,2,2) = oop_spouse(jage,1,2,1) ! no insurance
oop_spouse(jage,1,2,4) = oop_spouse(jage,1,2,1) ! no insurance
         end do   
close (159)


open(169, file='External_Parameters/OOP_spouse_old.txt')
         do jage=n_ret, n_age
read (169,*)  oop_spouse(jage,1,2,1) 
         end do 
close (169)


INC_TS(:)=0.0d0
open (232 , file='External_Parameters/parameters_inc_ts.txt')
read (232,*) INC_TS(1)
read (232,*) INC_TS(2)
read (232,*) INC_TS(3)
read (232,*) INC_TS(4)
close (232) 
 


PT_penalty(:)=0.0d0
stdze(:)=0.0d0
varzfix(:)=0.0d0
varzt(:)=0.0d0
p_zhc(:,:,:)=0.0d0 
  param_calib(:)=0.0d0
FT_fc(:,:)=0.d0   
eta(:,:)=0.0d0
beta_hat(:)=0.0d0
theta_beq=0.0d0
k_beq=0.0d0
sig_beq=0.0d0
cons_floor(:,:,:,:)=0.0d0
c_floor_adjust_mar(:)=0.0d0
c_floor_age_adj_young(:)=0.0d0
c_floor_age_adj_old(:)=0.0d0
In_ht(:,:)=0.0d0
delta_o(:,:)=0.0d0
delta=0.0d0
param_o_trans(:,:)=0.0d0
leis_mar(:,:,:)=0.0d0
gamma=0.0d0
zfix(:,:)=0.0d0
   
open(203, file='Internal_parameters/parameters_internal_preferences.txt')
read (203,*) alpha
read (203,*) sigma
read (203,*) cost_death
do je=1, n_e
read (203,*) beta_hat(je)
end do
 do je=1, n_e; do jh= 1, n_h
 read (203,*)  FT_fc( je, jh)   
 end do ; end do
 read (203,*) eta(2,1) 
 read (203,*) eta(2,2)
 read (203,*) eta(2,3) 
 read (203,*) eta(2,4)
 read (203,*) eta(2,5) 
do jo=1,5,2 ; do je=1,3      
  read (203,*)    leis_mar(je,jo,1) 
end do ; end do 
do jo=1,5,2 ; do je=1,3      
  read (203,*)    leis_mar(je,jo,2) 
end do ; end do 
read (203,*) theta_beq
read (203,*) k_beq
read (203,*) sig_beq
close (203)



if ((jExperimentNumber.eq.12).OR.(jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23)) THEN 
    eta(:,1)=eta(:,3)
    eta(:,2)=eta(:,3)
    eta(:,4)=eta(:,3)
end if 

    
    
    

open(204, file='Internal_parameters/parameters_internal.txt')
do je=1, n_e
read (204,*)	  PT_penalty(je)
read (204,*)      stdze(je)
read (204,*)      varzfix(je)
read (204,*)      varzt(je)
read (204,*)	  b0(je)	
read (204,*)	  b1(je)	
read (204,*)	  b2(je)
read (204,*)	  b3(je)	
read (204,*)	  b4(je)
read (204,*)	  b5(je)
end do
do je=1, n_e
    read (204,*)	  b6(je)
end do
do je=1, n_e
read (204,*) cons_floor(1,1,je, 1) 
read (204,*) cons_floor(1,1,je, 2) 
end do
do je=1, n_e
read (204,*) c_floor_adjust_mar(je) 
end do
do je=1, n_e
read (204,*) c_floor_age_adj_young(je)
read (204,*) c_floor_age_adj_old(je)
end do
do je=1, n_e ; do jfix=1, n_fix
read (204,*) In_ht(je,jfix)  
end do; end do
do je=1, n_e 
read (204,*) param_o_trans(je,1) 
read (204,*) param_o_trans(je,2) 
read (204,*) param_o_trans(je,3) 
read (204,*) param_o_trans(je,4) 
read (204,*) param_o_trans(je,5) 
read (204,*) param_o_trans(je,6) 
end do

 read (204,*) delta_o(1,1) 
  read (204,*) delta_o(2,1)
   read (204,*) delta_o(3,1)
    read (204,*) delta_o(1,2) 
     read (204,*) delta_o(2,2)
      read (204,*) delta_o(3,2)
read (204,*) zfix(1,1)
read (204,*) zfix(1,2) 
read (204,*) zfix(1,n_fix) 
read (204,*) zfix(2,1)
read (204,*) zfix(2,2)
read (204,*) zfix(2,n_fix) 
read (204,*) zfix(3,1)
read (204,*) zfix(3,2)
read (204,*) zfix(3,n_fix) 
close (204)



! Blacks 
zfix(1,1)=-0.48d0
zfix(1,2)=-0.07d0
zfix(1,3)=0.33d0
c_floor_adjust_mar(1) = c_floor_adjust_mar(1)/1.2d0 


do je=1,3 ; do jos=1,2
leis_mar(je,2,jos)=leis_mar(je,3,jos)
leis_mar(je,4,jos)=leis_mar(je,5,jos)
end do; end do

p_zhc(:,3,:) = p_zhc(:,2,:) 
p_zhc(:,4,:) = p_zhc(:,2,:)
p_zhc(:,5,:) = p_zhc(:,2,:)

do je=1, n_e
cons_floor(1:n_age,1,je, 1) =cons_floor(1,1,je, 1) 
cons_floor(1:n_age,1,je, 2) =cons_floor(1,1,je, 2) 
cons_floor(1:n_age,1,je, 3) =cons_floor(1,1,je, 2) 
end do

do je=1, n_e
do jage=1, n_ret-1
cons_floor(jage,2,je, 1) = cons_floor(jage,1,je, 1)* equiv(je,jage)  
cons_floor(jage,2,je, 2) = cons_floor(jage,1,je, 2)* equiv(je,jage)*(c_floor_adjust_mar(je)) 
cons_floor(jage,2,je, 3) = cons_floor(jage,1,je, 3)* equiv(je,jage)*(c_floor_adjust_mar(je)) 
    end do 
end do 

do je=1, n_e
do jage=n_ret, n_age
    do jh=1, n_h
cons_floor(jage,2,je, jh) = cons_floor(jage,1,je, jh)* equiv(je,jage)
end do
    end do 
end do 

do jage=21, n_ret-1; do je=1, n_e
   cons_floor(jage,2,je, 2:3) = cons_floor(jage,2,je, 2:3)+   cons_floor(jage,2,je, 2:3)* (dble(jage)-20.0d0)*c_floor_age_adj_old(je)
end do ; end do 

do jage=19, 1, -1; do je=1, n_e
   cons_floor(jage,2,je, 2:3) = cons_floor(jage,2,je, 2:3)-   cons_floor(jage,2,je, 2:3)* (20.0d0 - dble(jage))*c_floor_age_adj_young(je)
end do ; end do



if (jExperimentNumber.EQ.18) then 
    do je=1,n_e; do jage=1,40; do jh=2,3
    cons_floor(jage,1,je,jh) = cons_floor(jage,1,je,jh)*1.1d0
    end do; end do; end do
end if 
if (jExperimentNumber.EQ.19) then 
     do je=1,n_e;   do jage=1,40; do jh=2,3
    cons_floor(jage,1,je,jh) = cons_floor(jage,1,je,jh)*1.2d0
    end do; end do; end do
end if 


do je=1, n_e; do jh=1, n_h; do jmar=1,n_mar; do jage=1, n_age
cons_floor_at(jage,jmar,je,jh)    = cons_floor(jage,jmar,je,jh)*(1.0d0+taxc) 
end do; end do; end do; end do



hrs_offer(:)=0.0d0
hrs_offer(1)=0.0d0
hrs_offer(2)=hrs_pt 
hrs_offer(3)=hrs_pt 
hrs_offer(4)=hrs_ft 
hrs_offer(5)=hrs_ft 

wage_penalty(:,:)=0.0d0
wage_penalty(1,:) = 0.0d0
wage_penalty(2,1) = PT_penalty(1)
wage_penalty(2,2) = PT_penalty(2)
wage_penalty(2,3) = PT_penalty(3)
wage_penalty(3,1) = PT_penalty(1)
wage_penalty(3,2) = PT_penalty(2)
wage_penalty(3,3) = PT_penalty(3)
wage_penalty(4,:) = 1.0d0
wage_penalty(5,:) = 1.0d0

do je=1, n_e
b1(je)=b1(je)/2.d0
b2(je)=b2(je)/4.d0
b3(je)=b3(je)/8.d0
end do

zhc(:)=0.0d0

zhc(1)=1.0d0

HK(:)=0.0d0
do jx=0, n_hc 
HK(jx) = dble(jx)
end do

max_age_HC(:)=0.0d0
max_age_HC(1) = 0.0d0
max_age_HC(2) = zhc(n_zhc)  

do jage=3, n_ret-1
max_age_HC(jage) = (max_age_HC(jage-1) +1.0d0) *zhc(n_zhc)  
end do 

xfn(:, :)=0.0d0

do jage=1, n_ret-1
xfn(1, jage)=0.d0 
xfn(2, jage)=  max_age_HC(jage) * 0.22d0  
xfn(3, jage)=  max_age_HC(jage) * 0.45d0    
xfn(4, jage)=  max_age_HC(jage) * 0.65d0    
xfn(5, jage)= max_age_HC(jage)
end do

pen(:,:,:)=0.0d0 

open(187, file='Internal_parameters/SS_inc.txt')
do jfix=1, n_fix; do jhc_ret =1, n_hc_ret
read (187,*)	  pen(jhc_ret, 1, jfix)
end do; end do 
close(187)

hc_TS(:,:)=0.0d0
open(68, file='Internal_Parameters/parameters_hc_ts.txt')
do je =1, n_e
 read (68,*) hc_TS(1,je)
  read (68,*) hc_TS(2,je)
  end do
close (68)

p_fix(:,:)=1.0d0/3.0d0 

In_type1(:)=0.0d0
In_type1(1)=p_e(1)*(1.0d0/3.0d0)*In_ht(1,1)
In_type1(2)=p_e(1)*(1.0d0/3.0d0)*(1.0d0-In_ht(1,1))
In_type1(3)=p_e(1)*(1.0d0/3.0d0)*In_ht(1,2)
In_type1(4)=p_e(1)*(1.0d0/3.0d0)*(1.0d0-In_ht(1,2))
In_type1(5)=p_e(1)*(1.0d0/3.0d0)*In_ht(1,3)
In_type1(6)=p_e(1)*(1.0d0/3.0d0)*(1.0d0-In_ht(1,3))

In_type1(7)=p_e(2)*(1.0d0/3.0d0)*In_ht(2,1)
In_type1(8)=p_e(2)*(1.0d0/3.0d0)*(1.0d0-In_ht(2,1))
In_type1(9)=p_e(2)*(1.0d0/3.0d0)*In_ht(2,2)
In_type1(10)=p_e(2)*(1.0d0/3.0d0)*(1.0d0-In_ht(2,2))
In_type1(11)=p_e(2)*(1.0d0/3.0d0)*In_ht(2,3)
In_type1(12)=p_e(2)*(1.0d0/3.0d0)*(1.0d0-In_ht(2,3))

In_type1(13)=p_e(3)*(1.0d0/3.0d0)*In_ht(3,1)
In_type1(14)=p_e(3)*(1.0d0/3.0d0)*(1.0d0-In_ht(3,1))
In_type1(15)=p_e(3)*(1.0d0/3.0d0)*In_ht(3,2)
In_type1(16)=p_e(3)*(1.0d0/3.0d0)*(1.0d0-In_ht(3,2))
In_type1(17)=p_e(3)*(1.0d0/3.0d0)*In_ht(3,3)
In_type1(18)=p_e(3)*(1.0d0/3.0d0)*(1.0d0-In_ht(3,3))

Pr_o_25(:,:)=0.0d0
open(105, file='Internal_Parameters/parameters_offers_25.txt')
do je=1, n_e; do jo=1, n_o
   read (105,*) Pr_o_25(jo, je)  
end do; end do 
close (105)


! adjust these for blacks so slightly fewer start out with ESHI. 
Pr_o_25(2, 1) = 0.21d0
Pr_o_25(3, 1) = 0.03d0
Pr_o_25(4, 1) = 0.36d0
Pr_o_25(5, 1) = 0.40d0


if (jExperimentNumber.eq.22) THEN
param_o_trans(:,2) =param_o_trans(:,2)*2.0d0
param_o_trans(:,3) =param_o_trans(:,3)*2.0d0
param_o_trans(:,5) =param_o_trans(:,5)*2.0d0
param_o_trans(:,6) =param_o_trans(:,6)*2.0d0
end if 

Pr_o_i(:,:,:,:)=0.0d0
do je=1, n_e 
Pr_o_i(2,1,je,1) = (1.0d0 - param_o_trans(je,3)  ) * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(2,2,je,1) = (1.0d0 - param_o_trans(je,2)  ) * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(2,3,je,1) = (1.0d0 - param_o_trans(je,2)  ) * (1.0d0 - param_o_trans(je,4)  )
Pr_o_i(2,4,je,1) = (1.0d0 - param_o_trans(je,1)  ) * (1.0d0 - param_o_trans(je,6)  )
Pr_o_i(2,5,je,1) = (1.0d0 - param_o_trans(je,1)  ) * (1.0d0 - param_o_trans(je,4)  )

Pr_o_i(3,1,je,1) = (1.0d0 - param_o_trans(je,3)  ) *  param_o_trans(je,5)  
Pr_o_i(3,2,je,1) = (1.0d0 - param_o_trans(je,2)  ) *  param_o_trans(je,5)  
Pr_o_i(3,3,je,1) = (1.0d0 - param_o_trans(je,2)  ) *  param_o_trans(je,4)  
Pr_o_i(3,4,je,1) = (1.0d0 - param_o_trans(je,1)  ) *  param_o_trans(je,6)  
Pr_o_i(3,5,je,1) = (1.0d0 - param_o_trans(je,1)  ) *  param_o_trans(je,4)  

Pr_o_i(4,1,je,1) =  param_o_trans(je,3)   * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(4,2,je,1) =  param_o_trans(je,2)   * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(4,3,je,1) =  param_o_trans(je,2)   * (1.0d0 - param_o_trans(je,4)  )
Pr_o_i(4,4,je,1) =  param_o_trans(je,1)   * (1.0d0 - param_o_trans(je,6)  )
Pr_o_i(4,5,je,1) =  param_o_trans(je,1)   * (1.0d0 - param_o_trans(je,4)  )

Pr_o_i(5,1,je,1) =  param_o_trans(je,3)   *  param_o_trans(je,5)  
Pr_o_i(5,2,je,1) =  param_o_trans(je,2)   *  param_o_trans(je,5)  
Pr_o_i(5,3,je,1) =  param_o_trans(je,2)   *  param_o_trans(je,4)  
Pr_o_i(5,4,je,1) =  param_o_trans(je,1)   *  param_o_trans(je,6)  
Pr_o_i(5,5,je,1) =  param_o_trans(je,1)   *  param_o_trans(je,4)  

end do 

do jage=2,n_ret-1
Pr_o_i(:,:,:,jage) =   Pr_o_i(:,:,:,1)
end do

scale_param(:)=0.0d0
je=1
do jage=30, 35
    do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-29))*delta_o(je,1) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(1,jo,je,jage)=min(1.0d0,dble((jage-29))*delta_o(je,1) + Pr_o_i(1,jo,je,jage) ) 
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage)* scale_param(je))
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* scale_param(je))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
end do
end do 

do jage=36, n_ret-1
     do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-29))*delta_o(je,2) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(1,jo,je,jage)=min(1.0d0,dble((jage-29))*delta_o(je,2) + Pr_o_i(1,jo,je,jage) ) 
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage)* scale_param(je))
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* scale_param(je))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
end do
end do

scale_param(:)=0.0d0
je=2
do jage=23, 32 
    do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-22))*delta_o(je,1) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .85d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .85d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
    end do
end do 

do jage=33, n_ret-1
     do jo=1,5
scale_param(je)=  (1.0d0 - dble((jage-23))*delta_o(je,2) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .9d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .9d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
end do
end do

scale_param(:)=0.0d0
je=3
do jage=23, 32 
    do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-22))*delta_o(je,1) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))

Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .85d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .85d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
end do 
end do

do jage=33, n_ret-1 
     do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-23))*delta_o(je,2) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .9d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .9d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
end do
end do

prem_mcare(2) =prem_mcare(1) * 2.0d0 

prem_ephi(5,3)=prem_ephi(3,3) 

prem_ephi(3,2)=prem_ephi(3,1) 
prem_ephi(5,1)=prem_ephi(3,1) 
prem_ephi(5,2)=prem_ephi(3,1) 

prem_ephi(1,3)=prem_ephi(3,1)
prem_ephi(2,3)=prem_ephi(3,1)
prem_ephi(4,3)=prem_ephi(3,1)

if ((jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23)) THEN
    prem_ephi(1,1) = prem_ephi(3,1)
    prem_ephi(2,1) = prem_ephi(3,1)
    prem_ephi(4,1) = prem_ephi(3,1)
    prem_ephi(1,3)=prem_ephi(3,1)*2.0d0
    prem_ephi(2,3)=prem_ephi(3,1)*2.0d0
    prem_ephi(4,3)=prem_ephi(3,1)*2.0d0
    prem_ephi(1,2)=prem_ephi(3,1)
    prem_ephi(2,2)=prem_ephi(3,1)
    prem_ephi(3,2)=prem_ephi(3,1)
    prem_ephi(4,2)=prem_ephi(3,1)
    prem_ephi(5,2)=prem_ephi(3,1)
end if 

In_h(:, :,:)=0.0d0
In_r(:, :)=0.0d0

open(201, file='Internal_Parameters/parameters_H_25.txt')
do je=1, n_e; do jht=2,1,-1; do jh=1,n_h 
   read (201,*) In_h(je, jh,jht)
end do ; end do; end do
close (201)

open(207, file='External_Parameters/parameters_R25.txt')
   read (207,*) In_r(1, 2) 
   In_r(1, 1)=1.0d0 - In_r(1, 2) 
   In_r(2, :)=In_r(1, :) 
   In_r(3, :)=In_r(1, :)
close (207)

 param_h_cutoff_A(:)=0.0d0
 param_h_cutoff_B(:)=0.0d0
 param_h(:)=0.0d0
 neu_shocks(:,:,:)=1.0d0 
 neu_R=1.0d0
 Shock_type(:,:,:,:)=0.0d0
 treat_age_adj=0.0d0
 open(22, file='Internal_parameters/parameters_internal_health.txt')
  
   read (22,*) param_h_cutoff_A(1) 
   read (22,*) param_h_cutoff_A(2)
   read (22,*) param_h_cutoff_B(1) 
   read (22,*) param_h_cutoff_B(2)
read (22,*) param_h(1)   
read (22,*) param_h(2)
read (22,*) param_h(3) 
read (22,*) param_h(4)
read (22,*) param_h(5) 
read (22,*) param_h(6) 
read (22,*) param_h(7) 
read (22,*) param_h(8) 
read (22,*) param_h(9)
read (22,*) param_h(10)
read (22,*) param_h(11)
read (22,*) param_h(12)
read (22,*) param_h(13)
read (22,*) param_h(14) 
read (22,*) neu_shocks(1,1,1) 
read (22,*) neu_shocks(2,1,1) 
read (22,*) neu_shocks(3,1,1) 
read (22,*) neu_shocks(1,2,1) 
read (22,*) neu_shocks(2,2,1) 
read (22,*) neu_shocks(3,2,1) 
read (22,*) neu_shocks(1,3,1) 
read (22,*) neu_shocks(2,3,1) 
read (22,*) neu_shocks(3,3,1) 
read (22,*) neu_R(1) 
read (22,*) Shock_type(1,3,1,1) 
read (22,*) Shock_type(1,3,2,1) 
read (22,*) Shock_type(1,3,3,1) 
read (22,*) Shock_type(1,1,1,1) 
read (22,*) Shock_type(1,1,2,1) 
read (22,*) Shock_type(1,1,3,1) 
read (22,*) Shock_type(2,3,1,1) 
read (22,*) Shock_type(2,3,2,1) 
read (22,*) Shock_type(2,3,3,1) 
read (22,*) Shock_type(2,1,1,1) 
read (22,*) Shock_type(2,1,2,1) 
read (22,*) Shock_type(2,1,3,1) 
read (22,*) treat_age_adj 
close (22)

! Blacks 
neu_shocks(2,3,1) =0.5d0 
neu_shocks(3,3,1) =0.5d0
Shock_type(2,1,2,1)= 0.15d0 
Shock_type(2,1,3,1)=0.07d0
Shock_type(1,1,2,1)=0.49d0  
Shock_type(1,1,3,1)=.75d0




Shock_type(:,2,:,1) = Shock_type(:,1,:,1)
Shock_type(:,4,:,1) = Shock_type(:,1,:,1)
Shock_type(:,5,:,1) = Shock_type(:,3,:,1)
Shock_type(3,:,:,1) = 1.0d0 - Shock_type(1,:,:,1) -Shock_type(2,:,:,1)
do jage=1,20
Shock_type(:,:,:,jage) = Shock_type(:,:,:,1) 
end do 

do jage=21, n_ret-1
    Shock_type(:,:,:,jage) = Shock_type(:,:,:,1)
end do 

do jage=21, n_ret-1 ; do jh=2,3
Shock_type(1,1,jh,jage) = Shock_type(1,1,jh,1)- treat_age_adj*(dble(jage)-20.0d0) 
Shock_type(1,2,jh,jage) = Shock_type(1,2,jh,1)- treat_age_adj*(dble(jage)-20.0d0) 
Shock_type(1,4,jh,jage) = Shock_type(1,4,jh,1)- treat_age_adj*(dble(jage)-20.0d0) 

Shock_type(2,1,jh,jage) =( Shock_type(2,1,jh,jage) / (Shock_type(2,1,jh,jage) + Shock_type(3,1,jh,1))) * (1.0d0 - Shock_type(1,1,jh,jage) ) 
Shock_type(2,2,jh,jage) =( Shock_type(2,2,jh,jage) / (Shock_type(2,2,jh,jage) + Shock_type(3,2,jh,1))) * (1.0d0 - Shock_type(1,2,jh,jage) )
Shock_type(2,4,jh,jage) =( Shock_type(2,4,jh,jage) / (Shock_type(2,4,jh,jage) + Shock_type(3,4,jh,1))) * (1.0d0 - Shock_type(1,4,jh,jage) )
end do ; end do

Shock_type(3,:,:,:) = 1.0d0 - Shock_type(1,:,:,:) -Shock_type(2,:,:,:) 

if (jExperimentNumber.eq.7) THEN
Shock_type(2,:,:,:) =  1.0d0 - Shock_type(3,:,:,:)
Shock_type(1,:,:,:) = 0.0d0
END IF 

if (jExperimentNumber.eq.8) THEN
Shock_type(2,:,:,:) =  1.0d0 
Shock_type(1,:,:,:) = 0.0d0
Shock_type(3,:,:,:) = 0.0d0
END IF 

if (jExperimentNumber.eq.10) THEN
Shock_type(2,:,:,:) =  Shock_type(2,:,:,:) + Shock_type(3,:,:,:)
Shock_type(3,:,:,:) = 0.0d0
END IF 

if (jExperimentNumber.eq.30) THEN 
Shock_type(1,1,:,:) =  Shock_type(1,1,:,:) + Shock_type(3,1,:,:)
Shock_type(1,2,:,:) =  Shock_type(1,2,:,:) + Shock_type(3,2,:,:)
Shock_type(1,4,:,:) =  Shock_type(1,4,:,:) + Shock_type(3,4,:,:)

Shock_type(2,3,:,:) =  Shock_type(3,3,:,:)
Shock_type(2,5,:,:) =  Shock_type(3,5,:,:)

Shock_type(3,:,:,:) = 0.0d0
END IF 

if ((jExperimentNumber.eq.12).OR.(jExperimentNumber.eq.13).OR.(jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23)) THEN 
Shock_type(:,1,:,:) =  Shock_type(:,3,:,:)
Shock_type(:,2,:,:) =  Shock_type(:,3,:,:)
Shock_type(:,4,:,:) =  Shock_type(:,3,:,:)
END IF 

if ((jExperimentNumber.eq.9).OR.(jExperimentNumber.eq.17)) THEN
I_treat_Mcaid=1
else 
I_treat_Mcaid=0
end if


neu_R(2)=neu_R(1)
neu_R(4)=neu_R(1)

 x_h(:)=0.0d0  
 index_h=0.0d0
 
 Tp_h(:,:, :, :, :,:,:,:)=0.0d0
 trprob_h(:,:,:,:,:,:,:,:)=0.0d0
 
do jh = 1, 3 ; do je = 1, 3 ; do jage = 1,  n_ret-1 ; do jdp = 0, 1 ; do jdu = 0, 1 
      x_h = (/ 0.0d0, 0.0d0, DBLE(jdp), DBLE(jdu), 0.0d0, 0.0d0, 0.0d0, 0.0d0, DBLE(jage), DBLE((jage)**2), DBLE((jage)**3), 0.0d0, 0.0d0, 0.0d0 /)

    IF (jh==1) THEN
        x_h(1)=1.0d0
    ELSEIF (jh==2) THEN
        x_h(2)=1.0d0
    ENDIF
    IF (je==3) THEN
        x_h(12)=1.0d0
    ELSEIF (je==2) THEN
        x_h(13)=1.0d0
    ENDIF
    If (jage.ge.41) then
        x_h(14)=1.0d0
    end if 
    
    
        index_h=DOT_PRODUCT(param_h(:),x_h(:))

        trprob_h(jh,1,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1))) 
        trprob_h(jh,2,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        
        trprob_h(jh,1,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))
        trprob_h(jh,2,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))

end do; end do; end do; end do; end do

x_h(:) =0.0d0
index_h=0.0d0

do jh = 1, 3 ; do je = 1, 3 ; do jage = 1, n_ret-1 ; do jdp = 0, 1 ; do jdu = 0, 1 ; do jmc = 2, 5
    x_h = (/ 0.0d0, 0.0d0, DBLE(jdp), DBLE(jdu), 0.0d0, 0.0d0, 0.0d0, 0.0d0, DBLE(jage), DBLE((jage)**2), DBLE((jage)**3), 0.0d0, 0.0d0 ,0.0d0 /) 

    IF (jh==1) THEN
        x_h(1)=1.0d0
    ELSEIF (jh==2) THEN
        x_h(2)=1.0d0
    ENDIF
    IF (je==3) THEN
        x_h(12)=1.0d0
    ELSEIF (je==2) THEN
        x_h(13)=1.0d0
    ENDIF
    if (jmc==2) then
       x_h(5)=1.0d0 
       else if (jmc==3) then
       x_h(6)=1.0d0 
        else if (jmc==4) then
       x_h(7)=1.0d0 
        else if (jmc==5) then
       x_h(8)=1.0d0 
    end if 
    
        index_h=DOT_PRODUCT(param_h(1:13),x_h(1:13))

        trprob_h(jh,1,je,1,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        trprob_h(jh,2,je,1,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        
        trprob_h(jh,1,je,2,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))
        trprob_h(jh,2,je,2,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))

end do; end do; end do; end do; end do; end do

param_h_cutoff_A(1)=param_h_cutoff_A(1)+ 0.6d0
param_h_cutoff_B(1)=param_h_cutoff_B(1)+ 0.6d0

x_h(:) =0.0d0
index_h=0.0d0
 
do jh = 1, 3 ; do je = 1, 3 ; do jage = n_ret,  n_age ; do jdp = 0, 1 ; do jdu = 0, 1 
      x_h = (/ 0.0d0, 0.0d0, DBLE(jdp), DBLE(jdu), 0.0d0, 0.0d0, 0.0d0, 0.0d0, DBLE(jage), DBLE((jage)**2), DBLE((jage)**3), 0.0d0, 0.0d0, 0.0d0 /)

    IF (jh==1) THEN
        x_h(1)=1.0d0
    ELSEIF (jh==2) THEN
        x_h(2)=1.0d0
    ENDIF
    IF (je==3) THEN
        x_h(12)=1.0d0
    ELSEIF (je==2) THEN
        x_h(13)=1.0d0
    ENDIF
    If (jage.ge.41) then
        x_h(14)=1.0d0
    end if 
    
    
        index_h=DOT_PRODUCT(param_h(:),x_h(:))

        trprob_h(jh,1,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1))) 
        trprob_h(jh,2,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        
        trprob_h(jh,1,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))
        trprob_h(jh,2,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))

end do; end do; end do; end do; end do

trprob_h(:,3,:,:,:,:,:,:) = 1.0d0 - trprob_h(:,1,:,:,:,:,:,:) - trprob_h(:,2,:,:,:,:,:,:)

do jh_p = 1,n_h ; do jh = 1,n_h ; do je=1,n_e ; do jht=1,2 ;  do jage = 1, n_age; do jmc=1,4 
Tp_h(jh_p,jh,je,jht,jage,1,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,1)   
Tp_h(jh_p,jh,je,jht,jage,3,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,1)   

Tp_h(jh_p,jh,je,jht,jage,6,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,1) 
Tp_h(jh_p,jh,je,jht,jage,8,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,1) 

Tp_h(jh_p,jh,je,jht,jage,2,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,1) 
Tp_h(jh_p,jh,je,jht,jage,4,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,1) 

Tp_h(jh_p,jh,je,jht,jage,5,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,1) 
Tp_h(jh_p,jh,je,jht,jage,7,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,1) 

end do; end do; end do; end do; end do; end do

Tp_h(:,:,:,:,n_ret:n_age,:,1,:) = Tp_h(:,:,:,:,n_ret:n_age,:,2,:)

do jh_p = 1,n_h ; do jh = 1,n_h ;do je=1,n_e ; do jht=1,2 ;  do jage = 1, n_ret-1 ; do jmc=1,4 
Tp_h(jh_p,jh,je,jht,jage,1,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,1)   
Tp_h(jh_p,jh,je,jht,jage,3,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,jmc+1)   

Tp_h(jh_p,jh,je,jht,jage,6,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,jmc+1) 
Tp_h(jh_p,jh,je,jht,jage,8,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,jmc+1) 

Tp_h(jh_p,jh,je,jht,jage,2,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,jmc+1) 
Tp_h(jh_p,jh,je,jht,jage,4,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,jmc+1) 

Tp_h(jh_p,jh,je,jht,jage,5,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,jmc+1) 
Tp_h(jh_p,jh,je,jht,jage,7,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,jmc+1) 

end do; end do; end do; end do; end do; end do

shocks_risk(:,:, :, :, :)  =0.d0 

open(222, file='External_Parameters/parameters_health_shocks.txt')
do je=1,n_e; do jtype2=1, n_type2 ; do jh=1,n_h ; do jr=1,n_r;  do jage=1, n_age
   read (222,*) shocks_risk(je,jtype2,jh,jr,jage)
end do ; end do ; end do  ; end do ; end do
close (222) 

    if (jExperimentNumber.eq.2.OR.jExperimentNumber.eq.4) then

   shocks_risk(:,1, :, :, 1:n_ret-1) = shocks_risk(:,1, :, :, 1:n_ret-1) +    shocks_risk(:,5, :, :, 1:n_ret-1)
   shocks_risk(:,2, :, :, 1:n_ret-1) = shocks_risk(:,2, :, :, 1:n_ret-1) +    shocks_risk(:,6, :, :, 1:n_ret-1)
   shocks_risk(:,3, :, :, 1:n_ret-1) = shocks_risk(:,3, :, :, 1:n_ret-1) +    shocks_risk(:,7, :, :, 1:n_ret-1)
   shocks_risk(:,4, :, :, 1:n_ret-1) = shocks_risk(:,4, :, :, 1:n_ret-1) +    shocks_risk(:,8, :, :, 1:n_ret-1)
      
shocks_risk(:,5, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,6, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,7, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,8, :, :, 1:n_ret-1) = 0.0

    end if 
   
           if (jExperimentNumber.eq.3.OR.jExperimentNumber.eq.4) then
    shocks_risk(:,1, :, :, 1:n_ret-1) = shocks_risk(:,1, :, :, 1:n_ret-1) +    shocks_risk(:,2, :, :, 1:n_ret-1)
   shocks_risk(:,3, :, :, 1:n_ret-1) = shocks_risk(:,3, :, :, 1:n_ret-1) +    shocks_risk(:,4, :, :, 1:n_ret-1)
   shocks_risk(:,5, :, :, 1:n_ret-1) = shocks_risk(:,5, :, :, 1:n_ret-1) +    shocks_risk(:,6, :, :, 1:n_ret-1)
   shocks_risk(:,7, :, :, 1:n_ret-1) = shocks_risk(:,7, :, :, 1:n_ret-1) +    shocks_risk(:,8, :, :, 1:n_ret-1)
      
shocks_risk(:,2, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,4, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,6, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,8, :, :, 1:n_ret-1) = 0.0
           end if 

                  if (jExperimentNumber.eq.1.OR.jExperimentNumber.eq.4) then
   shocks_risk(:,1, :, :, 1:n_ret-1) = shocks_risk(:,1, :, :, 1:n_ret-1) +    shocks_risk(:,3, :, :, 1:n_ret-1)
   shocks_risk(:,2, :, :, 1:n_ret-1) = shocks_risk(:,2, :, :, 1:n_ret-1) +    shocks_risk(:,4, :, :, 1:n_ret-1)
   shocks_risk(:,5, :, :, 1:n_ret-1) = shocks_risk(:,5, :, :, 1:n_ret-1) +    shocks_risk(:,7, :, :, 1:n_ret-1)
   shocks_risk(:,6, :, :, 1:n_ret-1) = shocks_risk(:,6, :, :, 1:n_ret-1) +    shocks_risk(:,8, :, :, 1:n_ret-1)
      
shocks_risk(:,3, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,4, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,7, :, :, 1:n_ret-1) = 0.0
shocks_risk(:,8, :, :, 1:n_ret-1) = 0.0
           
       end if 

p_risk(:,:,:,:)=0.d0 
open(223, file='External_Parameters/parameters_health_R.txt')  
do jh=1,n_h ; do jr=1,n_r;  do jage=1, n_age
   read (223,*) p_risk(jr,2, jage,jh)
end do ; end do ; end do  
close (223) 

p_risk(:,1,:,:) = 1.0d0 - p_risk(:,2,:,:)  

 if (jExperimentNumber.eq.6) then 
p_risk(:,2, 1:n_ret-1,:)   = 0.0d0 
p_risk(:,1, 1:n_ret-1,:)   = 1.0d0 
In_r(:,2)=0.0d0 
In_r(:,1)=1.0d0 
end if 

oop(:, :, :,:,:)=0.0d0

open(198, file='External_Parameters/ME_ins_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (198,*) oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+1,1,5)=oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+2,1,5)=oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+3,1,5)=oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+4,1,5)=oop(jtype2, jh, jage,1,5)
end do ; end do ; end do  
close (198)   

oop(:, :, :,1,3) =oop(:, :,:,1,5) 

open(168, file='External_Parameters/ME_ins_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (168,*) oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+1,2,5)=oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+2,2,5)=oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+3,2,5)=oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+4,2,5)=oop(jtype2, jh, jage,2,5)
end do ; end do ; end do  
close (168)   

oop(:, :, :,2,3) =oop(:, :,:,2,5) 

open(45, file='External_Parameters/ME_no_ins_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (45,*) oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+1,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+2,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+3,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+4,1,1)=oop(jtype2, jh, jage,1,1)
end do ; end do ; end do  
close (45)   

oop(:, :, :,1,2) =oop(:, :,:,1,1) 
oop(:, :, :,1,4) =oop(:, :,:,1,1) 

open(46, file='External_Parameters/ME_no_ins_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (46,*) oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+1,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+2,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+3,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+4,2,1)=oop(jtype2, jh, jage,2,1)
end do ; end do ; end do  
close (46)   

oop(:, :, :,2,2) =oop(:, :,:,2,1) 
oop(:, :, :,2,4) =oop(:, :,:,2,1) 

if (jExperimentNumber == 11) then
    oop(:, :, 1:n_ret-1,:,:)=0.0d0
end if 

if (jExperimentNumber == 20) then 
    oop(:, :, 1:n_ret-1,:,:)=oop(:, :, 1:n_ret-1,:,:)*1.1d0
end if 
if (jExperimentNumber == 21) then 
    oop(:, :, 1:n_ret-1,:,:)=oop(:, :, 1:n_ret-1,:,:)*0.9d0
end if

if ((jExperimentNumber == 12).OR.(jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23).OR.(jExperimentNumber.eq.24)) then 
    oop(:, :, 1:n_ret-1,:,1)=oop(:, :, 1:n_ret-1,:,5)
    oop(:, :, 1:n_ret-1,:,2)=oop(:, :, 1:n_ret-1,:,5)
    oop(:, :, 1:n_ret-1,:,4)=oop(:, :, 1:n_ret-1,:,5)
end if 
    

if (jExperimentNumber == 23) then 
    oop(:, :, 1:n_ret-1,:,:)=oop(:, :, 1:n_ret-1,:,:)*1.05d0
end if

open(47, file='External_Parameters/ME_ret_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=41, 51,5
   read (47,*) oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+1,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+2,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+3,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+4,1,1)=oop(jtype2, jh, jage,1,1)
end do ; end do ; end do  
close (47)  

open(48, file='External_Parameters/ME_ret_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=41, 51,5
   read (48,*) oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+1,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+2,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+3,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+4,2,1)=oop(jtype2, jh, jage,2,1)
end do ; end do ; end do  
close (48) 

do jage=56,n_age,1
oop(:, :, jage,:,1)=oop(:, :, 55,:,1)
end do

oop(:, :, n_ret:n_age,:,2) = oop(:, :, n_ret:n_age,:,1)
oop(:, :, n_ret:n_age,:,3) = oop(:, :, n_ret:n_age,:,1)
oop(:, :, n_ret:n_age,:,4) = oop(:, :, n_ret:n_age,:,1)
oop(:, :, n_ret:n_age,:,5) = oop(:, :, n_ret:n_age,:,1)

open(49, file='External_Parameters/Charges_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (49,*) Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+1,1)=Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+2,1)=Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+3,1)=Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+4,1)=Charges(jtype2, jh, jage,1)
end do ; end do ; end do  
close (49) 

open(50, file='External_Parameters/Charges_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (50,*) Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+1,2)=Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+2,2)=Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+3,2)=Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+4,2)=Charges(jtype2, jh, jage,2)
end do ; end do ; end do  
close (50) 

open(67, file='External_Parameters/parameters_mc_ts.txt')
 read (67,*) MC_TS(1)
  read (67,*) MC_TS(2)
   read (67,*) MC_TS(3)
close (67)

Prob_MEcat(:)=0.0d0
Prob_MEcat(1) = 0.95d0   
Prob_MEcat(2) = 1.0d0- Prob_MEcat(1) 

phi(:,:,:)=0.0d0  

open(230, file='External_Parameters/sick_days.txt')
do jtype2=1, n_type2
read (230,*) phi(1,1,jtype2) 
phi(1,2,jtype2)=phi(1,1,jtype2)
end do 

do jtype2=1, n_type2
read (230,*) phi(1,3,jtype2) 
end do 

do jtype2=1, n_type2
read (230,*) phi(2,1,jtype2) 
phi(2,2,jtype2)=phi(2,1,jtype2) 
phi(3,1,jtype2)=phi(2,1,jtype2) 
phi(3,2,jtype2)=phi(2,2,jtype2) 
end do 

do jtype2=1, n_type2
read (230,*) phi(2,3,jtype2) 
phi(3,3,jtype2)=phi(2,3,jtype2) 
end do 

close (230)

phi_nw(:,:)=0.0d0
open(233, file='External_Parameters/sick_days_nw.txt')
do jh=2,n_h
do jtype2=1, n_type2
read (233,*) phi_nw(jh,jtype2) 
end do 
end do 
close (233)
phi_nw(1,:) =phi_nw(2,:) 

 if (jExperimentNumber.eq.5) then 
   phi(:,:,:)=0.0d0  
   phi_nw(:,:)=0.0d0
 end if 

leis(:,:,:,:)=0.0d0
do je=1, n_e ; do jh=1, n_h; do jtype2=1, n_type2
leis(je,1,jh,jtype2) = 1.0d0 - phi_nw(jh,jtype2) 
leis(je,2:3,jh,jtype2)=  1.0d0 - hrs_pt - .75*FT_fc( je, jh)  
leis(je,4:5,jh,jtype2)=  1.0d0 - hrs_ft - FT_fc( je, jh) 
end do ; end do ; end do

stdzt(:)=0.0d0
stdzt(1)=varzt(1)**0.5d0
stdzt(2)=varzt(2)**0.5d0
stdzt(3)=varzt(3)**0.5d0

Call Deterministic_z 

zt(:,:)=0.0d0
do je=1, n_e
zt(je,1)=-1.d0*stdzt(je) 
zt(je,n_zt)=1.d0*stdzt(je) 
Call discretize_normal(zt(je,1),zt(je,n_zt),n_zt,0.d0,stdzt(je),zt(je,:),p_zt(je,:,2)) 
Call discretize_normal(zt(je,1),zt(je,n_zt),n_zt,b6(je),stdzt(je),zt(je,:),p_zt(je,:,1)) 

end do

 z(:,:,:,:,:)=0.0d0
 base_earn(:,:,:,:,:,:)=0.0d0 
do je=1, n_e; do jfix=1, n_fix;  do jzt=1, n_zt 
do jx=0,n_hc;  do jh=1,n_h 
  z(jx,jh,je,jfix,jzt) = exp(zd(je,jx,jh)+zfix(je,jfix)+zt(je,jzt))  
  
  
 if (jExperimentNumber.eq.29) then
      z(jx,jh,je,jfix,jzt)= z(jx,jh,je,jfix,jzt)*1.02d0
  end if 
  
  base_earn(jx,jh,je,jfix,jzt,1)=0.0d0 
  base_earn(jx,jh,je,jfix,jzt,2:3)= z(jx,jh,je,jfix,jzt)*wage_penalty(2,je) 
  base_earn(jx,jh,je,jfix,jzt,4:5)= z(jx,jh,je,jfix,jzt)  
end do
end do; end do ; end do  ; end do

assets_initial(:)=0.0d0

as(:,:, :)=0.0d0 

as(1,:,:)=0.d0
as(2,:, :)=10000.d0/5200.d0
as(3,:, :)=20000.d0/5200.d0
as(4,:, :)=50000.d0/5200.d0
as(5,:, :)=75000.d0/5200.d0
as(6,:, :)=100000.d0/5200.d0 
as(7,:, :)=135000.d0/5200.d0  
as(8,:,:)=200000.d0/5200.d0 
as(9,:, :)=400000.d0/5200.d0 
as(10,:, :)=750000.d0/5200.d0
as(11,:, :)=1100000.d0/5200.d0
as(12,:, :)=1700000.d0/5200.d0

sav_ret(:)=0.0d0

open(231, file='External_Parameters/savings_grid_ret.txt')
do jsav=1, n_sav
   read (231,*) sav_ret(jsav)
end do 
close (231)
sav_ret(:)=sav_ret(:)/5200.0d0


sav(:,:)=0.0d0
open(252, file='External_Parameters/savings_grid.txt')
do jsav=1, n_sav
   read (252,*) sav(5,jsav)
end do 
close (252)
sav(5,:)=sav(5,:)/5200.0d0


open(253, file='External_Parameters/savings_grid_PT.txt')
do jsav=1, n_sav
   read (253,*) sav(1,jsav)
end do 
close (253)
sav(1,:)=sav(1,:)/5200.0d0
sav(2,:)=sav(1,:)
sav(3,:)=sav(1,:)
sav(4,:)=sav(1,:)

 ASP=0.0d0   
 AS_LB(:,:,:,:,:)=0.0d0
    do jas=1, n_as
        do jtype1=1, n_type1 
            do jage=1, n_ret-1
                do jo=1, n_o
                    do jsav=1,n_sav
       ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)  
        
        if (ASP.ge.as(n_as,jtype1,jage+1)) then 
            AS_LB(jas,jtype1,jage+1,jo,jsav)= n_as
        else
            jl_a=1
            ju_a=n_as
        do
            if (ju_a-jl_a <= 1) exit
            jm_a = int(( ju_a + jl_a ) / 2)
            if (ASP >= as(jm_a,jtype1,jage+1) ) then
	            jl_a = jm_a
            else
	            ju_a = jm_a
            end if
        end do
        end if 
        
       AS_LB(jas,jtype1,jage+1,jo,jsav)=jl_a  
                    end do
                end do
            end do 
        end do
    end do
    
 HK_temp=0.0d0
 HK_LB(:,:,:,:,:,:,:)=0.0d0
 
  
    do jage=1, n_ret-2 
    do jhc=1, n_hc_grid 
     do je=1, n_e
         do jtype2=1, n_type2
            do jh=1, n_h
                 jx_p_o(:)=0.0d0
                 jx_p_o(1) = xfn(jhc,jage)
                 jx_p_o(2) =min( (xfn(jhc,jage)+( max((hrs_pt - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft)))  ,  dble(n_hc) ) 
                 jx_p_o(3) = jx_p_o(2)
                 jx_p_o(4) = min( (xfn(jhc,jage)+( max((hrs_ft - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft)))  ,  dble(n_hc) )
                 jx_p_o(5) = jx_p_o(4)
        
        do jzhc=1,n_zhc 
            do jo=1, n_o 
          HK_temp  = max(0.0d0, min( jx_p_o(jo) *zhc(jzhc), dble(n_hc) )) 
         if (HK_temp.ge. xfn(n_hc_grid, jage+1) ) then 
            HK_LB(jhc,je,jtype2,jh,jzhc,jo,jage+1)= n_hc_grid  
        else
            jl_hc=1
            ju_hc=n_hc_grid
        do
            if (ju_hc-jl_hc <= 1) exit
            jm_hc = int(( ju_hc + jl_hc ) / 2)
            if (HK_temp >= xfn(jm_hc   , jage+1) ) then
	            jl_hc = jm_hc
            else
	            ju_hc = jm_hc
            end if
            HK_LB(jhc,je,jtype2,jh,jzhc,jo,jage+1)=jl_hc
        end do
        end if  
            end do     
        end do  

            end do 
         end do
     end do
    
    end do 
    end do
    
  prob_s_treat_new(:,:,:,:,:,:,:,:,:,:)=0.0d0  
  do jage=1,n_ret-1; do jo=1,n_o; do je=1,n_e; do jht=1, n_ht; do jh=1,n_h; do jr=1, n_r; do jtype2=1, n_type2 
       do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)
             call TYPE3_DET(jh,jr,jtype3)
                          
        
        prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) = Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(2,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,1)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(3,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,2)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(4,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,3)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(5,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,4)*p_risk(jr, jr_p,jage,jh) 
 
            
       end do; end do ; end do ; end do ; end do  ; end do ; end do
  end do; end do ; end do ; end do 
 
 prob_s_treat_new_64(:,:,:,:,:,:,:,:)=0.0d0    
    do je=1,n_e; do jht=1, n_ht; do jh=1,n_h; do jr=1, n_r; do jtype2=1, n_type2 
        do jh_p=1, n_h ; do jr_p=1, n_r ; do jtype2_p=1, n_type2 ; do jMEcat_p=1, n_MEcat 
                        call TYPE3_DET(jh_p,jr_p,jtype3_p)
             call TYPE3_DET(jh,jr,jtype3)
            
        
            prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,2,1)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(2,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,1)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(3,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,2)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(4,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,3)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(5,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,4)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)

        
        end do; end do ; end do ; end do ; end do  ; end do ; end do; end do ; end do
        
        
        
      prob_s_treat_new_ret(:,:,:,:,:,:,:,:,:,:)=0.0d0    
 do jage=n_ret, n_age-1;   do je=1,n_e; do jht=1, n_ht; do jh=1,n_h; do jr=1, n_r; do jtype2=1, n_type2 ; do jmar=1,n_mar
        do jh_p=1, n_h ; do jr_p=1, n_r ; do jtype2_p=1, n_type2 ; do jMEcat_p=1, n_MEcat ; do jmar_p=1,n_mar
                        call TYPE3_DET(jh_p,jr_p,jtype3_p)
             call TYPE3_DET(jh,jr,jtype3)   
prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) = Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) * Tp_mar(jmar_p,jmar,je,jage,1,1,1) 
        end do; end do ; end do ; end do ; end do  ; end do ; end do
        end do ; end do   ; end do ; end do ; end do
        
jage=n_age  
RET_PRIME(:,:,:, :,:,:,:,:,:)   =0.d0

do jhc_ret = 1, n_hc_ret 
 je = 1
do jfix=1,min(n_fix, I_state_fix)  

do jht=1,min(n_ht, I_state_ht)  
call TYPE1_DET(je,jht,jfix,jtype1)
do jas=1,n_as
do jMEcat=1, n_MEcat  
do jmar=1,n_mar
do jdu=1,n_du ; do jdp=1,n_dp ; do js=1,n_s 
call TYPE2_DET(jdu,jdp,js,jtype2)         
do jh=1, n_h ; do jr=1, n_r
    
Resources = interest*as(jas,jtype1,n_age) -prem_mcare(jmar) - oop(jtype2,jh,n_age,jMEcat,1) - oop_spouse(n_age,je,jmar,1) +pen(jhc_ret,  je,jfix) +  inc_spouse(n_age,je,1,jmar,jfix)  - max(0.0d0,  (( ((max(0.0d0, (interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret, je,jfix) +  inc_spouse(jage,je,1,jmar,jfix)-prem_mcare(jmar) - max(0.d0, oop(jtype2,jh,n_age,jMEcat,1) + oop_spouse(n_age,je,jmar,1)  - 0.075d0*((interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret,je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )))*5200.0d0) - tax_lambda(jmar) *( ((max(0.0d0, (interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret, je,jfix) +  inc_spouse(jage,je,1,jmar,jfix)-prem_mcare(jmar) - max(0.d0, oop(jtype2,jh,n_age,jMEcat,1) + oop_spouse(n_age,je,jmar,1)  - 0.075d0*((interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret,je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )))*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))
CON=0.0d0
ASP=0.0d0
TR=0.0d0
V_cf_trpay(:)=0.0d0
V_c_all_trpay(:) =0.0d0

if (Resources.le.cons_floor_at(jage, jmar,je,jh)) then
   CON= cons_floor(jage, jmar,je,jh)  
   call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
  RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=UTI + beta_hat(je)* ( cost_death + Bequest_val(0.0d0) )
else 
    
    call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
V_cf_trpay(1) = UTI + beta_hat(je)* ( cost_death + Bequest_val(Resources-cons_floor_at(jage, jmar,je,jh)) )
    call func_uti_ret(je,jh,jmar,jtype2,(Resources/(1.0d0+taxc)),UTI)
V_c_all_trpay(1) = UTI+ beta_hat(je)* ( cost_death + Bequest_val(0.0d0) )

    MAX_SAV=0.0d0
    INT_sav_upper_ret=0.0d0
    jsav_upper_ret=0.0d0
    MAX_SAV=Resources-cons_floor_at(jage, jmar,je,jh) -as(jas,jtype1,n_age)
    
                            jsav_upper_ret=0
          do while ( MAX_SAV .ge.  sav_ret(jsav_upper_ret+1))
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav) exit
          end do 
          
        maxB=V_cf_trpay(1)
    do jsav= jsav_upper_ret , 1, -1
    CON= ( Resources -  as(jas,jtype1,n_age) - sav_ret(jsav) )/(1.0d0+taxc)
    ASP=  as(jas,jtype1,n_age)+  sav_ret(jsav)
    if (ASP.lt.0.d0)  exit
    if (CON.lt.cons_floor(jage, jmar,je,jh)) then 
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
    end if  
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
	    RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=    UTI+ beta_hat(je)* ( cost_death + Bequest_val(ASP) )
  
        
        if (RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar).lt.maxB) exit
			  	    maxB=RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)                  
   
    end do 
    
     RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=maxB 
     
    if (RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar).eq.0.0d0) then
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
    end if 
    
  
        RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=max(RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar),  V_c_all_trpay(1)) 
   
    
    if (RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar).eq.0.0d0) then
            PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
        end if 
        
end if 

                      
end do; end do 
end do 
end do; end do; end do 
end do 
end do 
end do 
end do 
end do 


do jage=n_age-1,n_ret,-1  
   PRINT *,'Current Age:',jage
   RET(:,:,:,:,:,:,:,:)=0.d0   
jfix=1 
je =1  
do jht=1, min(n_ht, I_state_ht)    
    
call TYPE1_DET(je,jht,jfix,jtype1)

do jdu=1,n_du ; do jdp=1,n_dp 
    js=1 
    
call TYPE2_DET(jdu,jdp,js,jtype2) 
do jh=1, n_h ; do jr=1, n_r  
    
    call TYPE3_DET(jh,jr,jtype3)
    
do js=1,n_s 
    call TYPE2_DET(jdu,jdp,js,jtype2) 
do jas=1,n_as
do jMEcat=1, n_MEcat  
do jmar=1,n_mar 

do jhc_ret = 1, n_hc_ret 

   TR=0.d0
   CON=0.d0
   ASP=0.d0
   UTI=0.d0
   TR_indic=0
   MAX_SAV=0.0d0
   V_cf_trpay(:)=0.0d0
   V_c_all_trpay(:) =0.0d0
    
income = max(0.0d0, (interest-1.0d0)*as(jas,jtype1,jage) + pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -  prem_mcare(jmar)  - max(0.d0,  oop(jtype2, jh, jage,jMEcat,1) +  oop_spouse(jage,je,jmar,1) - 0.075d0*((interest-1.0d0)*as(jas,jtype1,jage) + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )) 
base_2 = prem_mcare(jmar) + oop(jtype2, jh, jage,jMEcat,1) + oop_spouse(jage,je,jmar,1)  + max(0.0d0,  (( (income*5200.0d0) - tax_lambda(jmar) *( (income*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) 

if (((interest*as(jas,jtype1,jage)+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) .lt. (cons_floor_at(jage, jmar,je,jh)) ) then                                                                                                                                      
    TR= cons_floor_at(jage, jmar,je,jh) - (interest*as(jas,jtype1,jage)+ pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2 ) 
    CON=cons_floor(jage, jmar,je,jh)
    TR_indic=1 
    MAX_SAV=0.0d0 
else
   TR=0.d0
   CON=0.d0
   TR_indic=0
   MAX_SAV=  interest*as(jas,jtype1,jage)  + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) - base_2 -cons_floor_at(jage, jmar,je,jh) -as(jas,jtype1,jage)
end if 

if (TR_indic.eq.1) then
    RET_PRIME_INT(:,:,:,:,:)=0.0d0
    EXPECTED_RET_PRIME=0.d0 
    V_death=0.0d0
    V_death=cost_death +Bequest_val(0.0d0) 
    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
        call TYPE3_DET(jh_p,jr_p,jtype3_p)
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) *( p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  + (1.0d0- p_surv(jh_p,je,jmar_p, jage+1))  * V_death )
    end do ; end do  ; end do; end do; end do
 
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
        RET(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,jmar)=UTI+ beta_hat(je)*EXPECTED_RET_PRIME 
    
else if (TR_indic.eq.0) then
    
    RET_PRIME_INT(:,:,:,:,:)=0.0d0
    EXPECTED_RET_PRIME=0.d0 
    ASP=0.0d0
    ASP=MAX_SAV+as(jas,jtype1,jage)
    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
        call TYPE3_DET(jh_p,jr_p,jtype3_p)
                call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) *( p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)  + (1.0d0- p_surv(jh_p,je,jmar_p, jage+1))  * (cost_death +Bequest_val(ASP) ) )
    end do ; end do  ; end do; end do; end do
    
    call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
  V_cf_trpay(1) = UTI + beta_hat(je)* EXPECTED_RET_PRIME

  EXPECTED_RET_PRIME=0.0d0
      do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
          call TYPE3_DET(jh_p,jr_p,jtype3_p)
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) *( p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  + (1.0d0- p_surv(jh_p,je,jmar_p, jage+1))  * (cost_death +Bequest_val(0.0d0) ) )
    end do ; end do  ; end do; end do; end do
    call func_uti_ret(je,jh,jmar,jtype2,((((interest*as(jas,jtype1,jage)+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) )/(1.0d0+taxc)),UTI) 
V_c_all_trpay(1) = UTI+ beta_hat(je)* EXPECTED_RET_PRIME

                     jsav_upper_ret=0
          do while ( MAX_SAV .ge.  sav_ret(jsav_upper_ret+1))
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav) exit
          end do 
          
        maxB=V_cf_trpay(1) 
    do jsav= jsav_upper_ret , 1, -1
    TR=0.d0
    CON= ( ((interest-1.0d0))*as(jas,jtype1,jage) + inc_spouse(jage,je,1,jmar,jfix) + pen(jhc_ret,  je,jfix) -base_2  - sav_ret(jsav) )/(1.0d0+taxc)
    ASP=  as(jas,jtype1,jage)+  sav_ret(jsav)
    
    if (ASP.lt.0.d0)  exit
    if (CON.lt.cons_floor(jage, jmar,je,jh)) then 
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
    end if  
           
            EXPECTED_RET_PRIME=0.d0
            RET_PRIME_INT(:,:,:,:,:)=0.0d0
             do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
                 call TYPE3_DET(jh_p,jr_p,jtype3_p)
                call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )   
            EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) * (    p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)  +   (1.0d0-p_surv(jh_p,je,jmar_p, jage+1))* (cost_death + Bequest_val(ASP) ))
             end do; end do; end do  ; end do  ; end do 
             call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
	        RET(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,jmar)=UTI+ beta_hat(je)*EXPECTED_RET_PRIME    
       
            if (RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar).lt.maxB) exit
			  	    maxB=RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)                  
		        
                
    end do 
     RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)=maxB 
    
    if (RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar).eq.0.0d0) then
        RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)=V_c_all_trpay(1)   
    else 
        RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)=max(RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar), V_c_all_trpay(1)) 
    end if 
 
end if  
   	       
RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,jage,jmar)=RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)

end do 
end do 
end do 
end do 
end do; end do; end do; end do ;  end do  	 
end do 

end do 

RET_PRIME(:,:,:, 3,:,:,:,:,:)=RET_PRIME(:,:,:, 1,:,:,:,:,:)
RET_PRIME(:,:,:, 5,:,:,:,:,:)=RET_PRIME(:,:,:, 1,:,:,:,:,:)

RET_PRIME(:,:,:, 4,:,:,:,:,:)=RET_PRIME(:,:,:, 2,:,:,:,:,:)
RET_PRIME(:,:,:, 6,:,:,:,:,:)=RET_PRIME(:,:,:, 2,:,:,:,:,:)

RET_PRIME(:,:,:, 9,:,:,:,:,:)=RET_PRIME(:,:,:, 7,:,:,:,:,:)
RET_PRIME(:,:,:, 11,:,:,:,:,:)=RET_PRIME(:,:,:, 7,:,:,:,:,:)

RET_PRIME(:,:,:, 10,:,:,:,:,:)=RET_PRIME(:,:,:, 8,:,:,:,:,:)
RET_PRIME(:,:,:, 12,:,:,:,:,:)=RET_PRIME(:,:,:, 8,:,:,:,:,:)

RET_PRIME(:,:,:, 15,:,:,:,:,:)=RET_PRIME(:,:,:, 13,:,:,:,:,:)
RET_PRIME(:,:,:, 17,:,:,:,:,:)=RET_PRIME(:,:,:, 13,:,:,:,:,:)

RET_PRIME(:,:,:, 16,:,:,:,:,:)=RET_PRIME(:,:,:, 14,:,:,:,:,:)
RET_PRIME(:,:,:, 18,:,:,:,:,:)=RET_PRIME(:,:,:, 14,:,:,:,:,:)


!************************************************************************************************************************************
! AGE 64
!************************************************************************************************************************************

earnings_interpol=0.0d0 
EXP_L_prime(:,:,:,:,:,:,:,:)=0.0d0
B(:,:,:,:,:,:,:,:)=0.0d0

jage=n_ret-1   
PRINT *,'Current Age:',jage     

je =1
do jfix=1, min(n_fix , I_state_fix)
do jht=1, min(n_ht, I_state_ht)
    call TYPE1_DET(je,jht,jfix,jtype1)
    
EXP_L_prime_fix(:,:,:,:,:,:,:)=0.0d0
B_fix(:,:,:,:,:,:,:,:,:,:)=0.0d0
B_1(:,:,:,:,:,:,:,:,:)=0.0d0
B_2(:,:,:,:,:,:,:,:,:)=0.0d0
EXP_L_prime_simple(:,:,:,:,:,:,:)=0.0d0
    
do jhc =1, n_hc_grid     
  
do jtype4=1,n_type4
                if (jtype4.eq.1) then
                jmar=1
                jos=1
            else if (jtype4.eq.2) then
                jmar=2
                jos=1
            else if (jtype4.eq.3) then
                jmar=2
                jos=2
            end if 
do jh=1, n_h ; do jr=1, n_r  
    call TYPE3_DET(jh,jr,jtype3)
do jas=1,n_as ; do jzt=1, n_zt    
do jdu=1,n_du ; do jdp=1,n_dp ; do js=1,n_s
    call TYPE2_DET(jdu,jdp,js,jtype2) 

do jMEcat=1, n_MEcat  
     
do jo=1, n_o
    
         if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
        
  actual_earn =  0.0d0
  income_bt=0.0d0
  tax_SS=0.0d0
  tax_med=0.0d0
  income_pay(:)=0.0d0
  Total_tax_pay(:)=0.0d0
  base_2_pay(:)=0.0d0
  base_1_pay(:)=0.0d0
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay(:)=0.0d0
  V_c_all_trpay(:)=0.0d0
             SAVE_RET_PRIME_INT(:,:,:,:,:,:)=0.0d0
           SAVE_EXPECTED_B_PRIME_trpay(:,:)=0.0d0
  TR_indic=0
  
    if (jo.eq.1)  then 
      earnings_interpol=0.0d0
  else 
      call InterpolateScalar ( HK(0 : n_hc), base_earn(0:n_hc,jh,je,jfix,jzt,jo) , xfn(jhc,jage)  , earnings_interpol ) 
  end if 

earnings_capacity=0.0d0

 if (xfn(jhc,jage).le.hc_TS(1,je)) then 
    jhc_ret=1
else if (xfn(jhc,jage).le.hc_TS(2,je)) then 
    jhc_ret=2
else 
    jhc_ret=3
end if  
  
actual_earn = max(0.0d0 , ( earnings_interpol * (hrs_offer(jo) -phi(je,jh,jtype2)  )) ) 
income_bt= ((interest-1.0d0))*as(jas,jtype1,jage) +  actual_earn  

if ((jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23)) THEN 
 
if (jo.ne.1) then 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)))
end if 
    
else 
    

if ((jo.eq.3).OR.(jo.eq.5)) then 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)))
end if 

end if 
     
if (income_bt.le.INC_TS(1)) then 
jinc=1
else if (income_bt.le.INC_TS(2)) then 
jinc=2
else if (income_bt.le.INC_TS(3)) then 
jinc=3
else if (income_bt.le.INC_TS(4)) then 
jinc=4
else 
jinc=5
end if


income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo) + oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  
  base_2_pay(1) =  prem_ephi(jo,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo) +  oop_spouse(jage,je,jtype4,jo) +Total_tax_pay(1)  
  base_1_pay(1) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        

        EXPECTED_B_PRIME_C_FLOOR=0.d0  
        RET_PRIME_INT(:,:,:,:,:)=0.0d0
        V_death=0.0d0
        V_death= cost_death + Bequest_val(0.0d0) 

        do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r 
        do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar    
            call TYPE3_DET(jh_p,jr_p,jtype3_p)
        EXPECTED_B_PRIME_C_FLOOR = EXPECTED_B_PRIME_C_FLOOR + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) * Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
        end do; end do                                                                                 
        end do; end do ; end do
        
     EXPECTED_B_PRIME_C_FLOOR_NT =0.0d0
        do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r 
        do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar   
              call TYPE3_DET(jh_p,jr_p,jtype3_p)
        EXPECTED_B_PRIME_C_FLOOR_NT = EXPECTED_B_PRIME_C_FLOOR_NT + prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) * Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
        end do; end do                                                                                 
        end do; end do ; end do
    
 if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then  
    CON=cons_floor(jage, jmar,je,jh)
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON,UTI)   
    B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = UTI+ beta_hat(je)*EXPECTED_B_PRIME_C_FLOOR
  
 else 
    MAX_SAV_pay(1) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  -cons_floor_at(jage, jmar,je,jh) 
    
              jsav_upper_pay(1)=0
          do while ( MAX_SAV_pay(1) .ge.  sav(jo,jsav_upper_pay(1)+1))
           jsav_upper_pay(1)=jsav_upper_pay(1) +1
           if (jsav_upper_pay(1).EQ.n_sav) exit
          end do 
          
               if (jsav_upper_pay(1).lt.1) then
                 PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
               end if 
               
            
      
     V_cf_trpay(1)=0.0d0
     ASP= MAX_SAV_pay(1) +as(jas,jtype1,jage)
     V_death=0.0d0
     V_death=cost_death + Bequest_val(ASP) 
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)      
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
 
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))
          
          
          do while ((MAX_SAV_pay(1) +as(jas,jtype1,jage)) .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      
      
do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
    call TYPE3_DET(jh_p,jr_p,jtype3_p)
         if (ASP.eq.0.0d0) then
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else if (( m ).eq.n_as) then 
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                slope_as=0.0d0
                inter_as=0.0d0
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(m+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((m+1),jtype1,jage+1) -  as((m),jtype1,jage+1)  )
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((m),jtype1,jage+1)
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
    
    
EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
end do; end do  ; end do;   end do; end do
  V_cf_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(1)

     V_c_all_trpay(1) =0.0d0
     call func_uti(jage,1,je,jh,jo,jtype4,jtype2, (base_1_pay(1)/(1.0d0+taxc)), UTI)  
     V_c_all_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR 

     

           maxB_trpay(1)=V_cf_trpay(1) 
do jsav= jsav_upper_pay(1) , 1, -1 
                 
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
         if (ASP.le.0.d0)  exit 
    V_death=0.0d0
    V_death= cost_death + Bequest_val(ASP) 
    CON_pay(1)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav(jo,jsav) )/(1.0d0+taxc)
    
         if (CON_pay(1).lt.cons_floor(jage, jmar,je,jh)) then 
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
         end if 
    
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON_pay(1),UTI)
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
    call TYPE3_DET(jh_p,jr_p,jtype3_p)
         if (( AS_LB(jas,jtype1,jage+1,jo,jsav) ).eq.n_as) then 
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                slope_as=0.0d0
                inter_as=0.0d0
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jage+1) -  as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)  )
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
        SAVE_RET_PRIME_INT(jsav,jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)
        
EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
end do; end do  ; end do;   end do; end do

SAVE_EXPECTED_B_PRIME_trpay(jsav,1)=EXPECTED_B_PRIME_trpay(1) 
     
B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI + beta_hat(je)*EXPECTED_B_PRIME_trpay(1)
           
        if (B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).lt.maxB_trpay(1)) exit
              maxB_trpay(1)=B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)   
                                   
end do   
        B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(1),  V_c_all_trpay(1)) 
   
end if   
 
if (jtype2.eq.1) then      
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
  
else 
  
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  
  base_2_pay(2) =  prem_ephi(jo,jtype4) +  oop_spouse(jage,je,jtype4,jo) + Total_tax_pay(2)  
  base_1_pay(2) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        
   
 if ( base_1_pay(2) .gt. (cons_floor_at(jage, jmar,je,jh))) then 
 
   MAX_SAV_pay(2) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(2)  -cons_floor_at(jage, jmar,je,jh) 
   
   
             jsav_upper_pay(2)=0
          do while ( MAX_SAV_pay(2) .ge.  sav(jo,jsav_upper_pay(2)+1))
           jsav_upper_pay(2)=jsav_upper_pay(2) +1
           if (jsav_upper_pay(2).EQ.n_sav) exit
          end do 
          
     V_cf_trpay(2:3)=0.0d0
     ASP= MAX_SAV_pay(2) +as(jas,jtype1,jage)  
     V_death=0.0d0
     V_death=cost_death + Bequest_val(ASP) 
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
 
       if (ASP .lt. as(AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2)),jtype1,jage+1)) then 
          pause 
      else  
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2))
          
          do while (ASP .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      end if 
      
do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
    call TYPE3_DET(jh_p,jr_p,jtype3_p)

         if (ASP.eq.0.0d0) then
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else if (( m ).eq.n_as) then 
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                slope_as=0.0d0
                inter_as=0.0d0
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(m+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((m+1),jtype1,jage+1) -  as((m),jtype1,jage+1)  )
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((m),jtype1,jage+1)
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
    
EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) +  prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)

EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) +  prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo)*(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)

end do; end do  ; end do;   end do; end do
   call func_uti(jage,2,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
   V_cf_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(2)
  call func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
   V_cf_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(3)

     V_c_all_trpay(2) =0.0d0
     call func_uti(jage,2,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI)  
     V_c_all_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR 
     V_c_all_trpay(3) =0.0d0
     V_death=0.0d0
     V_death= cost_death + Bequest_val(0.0d0) 
     call func_uti(jage,3,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI) 
     V_c_all_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 

      
     maxB_trpay(2)=V_cf_trpay(2) 
     maxB_trpay(3)=V_cf_trpay(3)
do jsav= jsav_upper_pay(2) , 1, -1 
                    
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
    if (ASP.le.0.d0)  exit 
    
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
    
    CON_pay(2)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(2) -sav(jo,jsav) )/(1.0d0+taxc)
        
        
         if (CON_pay(2).lt.cons_floor(jage, jmar,je,jh)) then 
          PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
          pause
         end if 

                    call  func_uti(jage,2,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(2) = UTI
                    call  func_uti(jage,3,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(3) = UTI     
                                      
if ((SAVE_EXPECTED_B_PRIME_trpay(jsav,1)).ne.(0.0d0)) then 
    B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(2) + beta_hat(je)*SAVE_EXPECTED_B_PRIME_trpay(jsav,1)
  
    EXPECTED_B_PRIME_trpay(3)=0.0d0
    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
        call TYPE3_DET(jh_p,jr_p,jtype3_p)
    EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) + prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* SAVE_RET_PRIME_INT(jsav,jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
    end do; end do  ; end do;   end do; end do
    B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
else 
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
     do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
         call TYPE3_DET(jh_p,jr_p,jtype3_p)
         if (( AS_LB(jas,jtype1,jage+1,jo,jsav) ).eq.n_as) then 
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                slope_as=0.0d0
                inter_as=0.0d0
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jage+1) -  as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)  )
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    

            EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death )
            EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) + prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death )
     
     end do; end do  ; end do;   end do; end do
     
           B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(2) + beta_hat(je)*EXPECTED_B_PRIME_trpay(2)
           B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
end if            
 

     if ((B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).lt.maxB_trpay(2)).AND.(B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).LT.maxB_trpay(3)))  exit  
         
              if (B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).ge.maxB_trpay(2)) then                              
              maxB_trpay(2)=B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)    
              end if   
              
              if (B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).ge.maxB_trpay(3)) then                              
              maxB_trpay(3)=B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  
              end if 
              
         
end do   

        B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max( maxB_trpay(2),  V_c_all_trpay(2)) 
        B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max( maxB_trpay(3),  V_c_all_trpay(3)) 

    
        
               
 else 
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
    call  func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  =   UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 
 
 if ((I_treat_Mcaid.EQ.1).AND.( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) ) THEN
      B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = max( B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) , B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)) 
 END IF 
 
 
 
 end if   
end if  

if (jtype2.eq.1) then
B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) )
else 
    if ((Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1)).NE.0.0d0) then 
B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) = (Shock_type(2,jo,jh,n_ret-1)/(Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1)))* (max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) ))    + (Shock_type(1,jo,jh,n_ret-1)/(Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1)))* B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1)

    end if 
end if

   
if (jExperimentNumber.eq.13) THEN 
    B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) =  B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1)

    else 
B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1),B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) ,B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) )
    end if 
    
end do 
end do 
end do; end do; end do 

        EXPECTED_B(:)=0.0d0 
    do jo=1, n_o
            do jtype2=1, n_type2 
            do jMEcat=1, n_MEcat
                       
                prob_s=0.0d0
                prob_s= shocks_risk(je,jtype2,jh,jr, jage)*Prob_MEcat(jMEcat) 
   EXPECTED_B(jo)= EXPECTED_B(jo)+  prob_s * ((Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1))* B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1)     + Shock_type(3,jo,jh,n_ret-1) *B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) )
            end do 
            end do  
            B(jo,jtype4,jzt,jas,jtype3,jtype1,jhc,n_ret-1) = EXPECTED_B(jo)
    end do  
  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,:,jtype4)=EXPECTED_B(1) ! initialize to the non-employed value
    
    If   (EXPECTED_B(5).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,5,jtype4)=EXPECTED_B(5)
    end if
  
    If   (EXPECTED_B(4).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,4,jtype4)=EXPECTED_B(4)
    end if  
    
    If   (EXPECTED_B(3).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,3,jtype4)=EXPECTED_B(3)
    end if
  
    If   (EXPECTED_B(2).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,2,jtype4)=EXPECTED_B(2)
    end if 

end do; end do; end do ; end do   
end do 
end do 


do jo=1,n_o  ; do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid  
        if ((jtype3.eq.1).OR.(jtype3.eq.2)) then 
        jh=1
        else if ((jtype3.eq.3).OR.(jtype3.eq.4)) then
    jh=2
    else if ((jtype3.eq.5).OR.(jtype3.eq.6)) then
        jh=3
    end if 
    
   EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,:,1) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,1) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,1)  + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,1) * p_zt(je,jzt,1)
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,1) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,1)  + p_zt(je,jzt,1) * (EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,2) *prob_os(n_ret-1,je,jh,1,jfix) + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,3) *prob_os(n_ret-1,je,jh,2,jfix))
    end do
   
       EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,:,2) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,2) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,2)  + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,1) * p_zt(je,jzt,2)
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,2) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,2)  + p_zt(je,jzt,2) * (EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,2) *prob_os(n_ret-1,je,jh,1,jfix) + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,3) *prob_os(n_ret-1,je,jh,2,jfix))
    end do
    
end do; end do; end do; end do



do jo=1,n_o  ; do jmar =1,n_mar ; do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid 
EXP_L_prime(jas,jtype1,jtype3,jhc,n_ret-1,jo,jmar,1) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,jmar,1)
EXP_L_prime(jas,jtype1,jtype3,jhc,n_ret-1,jo,jmar,2) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,jmar,2)
end do; end do; end do; end do; end do


end do; end do 

 
!******************************************************************************************************************************************
! AGES 25-63  
!******************************************************************************************************************************************

je = 1
do jfix=1, min(n_fix , I_state_fix)
do jht=1, min(n_ht, I_state_ht)
call TYPE1_DET(je,jht,jfix,jtype1) 
B_fix(:,:,:,:,:,:,:,:,:,:)=0.0d0
EXP_L_prime_fix(:,:,:,:,:,:,:)=0.0d0

B_1(:,:,:,:,:,:,:,:,:)=0.0d0
B_2(:,:,:,:,:,:,:,:,:)=0.0d0

ageloop: do jage=n_ret-2, 1, -1   
    
      
  if ((jExperimentNumber.eq.25).AND.(jage.eq.6)) then
       base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)*1.02d0
  else if ((jExperimentNumber.eq.25).AND.(jage.eq.5)) then
       base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)/1.02d0
  end if 
  
  if ((jExperimentNumber.eq.26).AND.(jage.eq.16)) then
      base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)*1.02d0
  else if ((jExperimentNumber.eq.26).AND.(jage.eq.15)) then  
       base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)/1.02d0
  end if 
  
    if ((jExperimentNumber.eq.27).AND.(jage.eq.26)) then
      base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)*1.02d0
  else if ((jExperimentNumber.eq.27).AND.(jage.eq.25)) then  
       base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)/1.02d0
  end if 
  
      if ((jExperimentNumber.eq.28).AND.(jage.eq.36)) then
      base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)*1.02d0
  else if ((jExperimentNumber.eq.28).AND.(jage.eq.35)) then  
       base_earn(:,:,:,:,:,:)= base_earn(:,:,:,:,:,:)/1.02d0
  end if 
  
  
hkloop: do jhc =1, n_hc_grid     
    do jh=1, n_h 
do jdu=1,n_du ; do jdp=1,n_dp ; do js=1,n_s 
    call TYPE2_DET(jdu,jdp,js,jtype2)
    


jx_p_o(:)=0.0d0
jx_p_o(1) = xfn(jhc,jage)
jx_p_o(2) = min((xfn(jhc,jage)+( max((hrs_pt - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft))) ,  dble(n_hc) )
jx_p_o(3) = jx_p_o(2)
jx_p_o(4) = min( (xfn(jhc,jage)+( max((hrs_ft - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft))) ,  dble(n_hc) )
jx_p_o(5) = jx_p_o(4)

do jzhc=1,n_zhc  
  HK_INTERPOL(jzhc,:)  =  max(0.0d0, min( jx_p_o(:) *zhc(jzhc), dble(n_hc) )) 
end do  

     
do jr=1, n_r 
    call TYPE3_DET(jh,jr,jtype3)  ! Gives Type3 given jh and jr
do jas=1,n_as  
do jtype4=1,n_type4 
!find jmar and jos    
            if (jtype4.eq.1) then
                jmar=1
                jos=1
            else if (jtype4.eq.2) then
                jmar=2
                jos=1
            else if (jtype4.eq.3) then
                jmar=2
                jos=2
            end if    
do jMEcat=1, n_MEcat  
  do jo=1,n_o 
      
      
      ! jo_lagged is 1 if not working, 2 if working. this is used when calculating the expected value tomorrow which depends on the current working status since wages tomorrow depend on lagged employmend. 
      if (jo.eq.1) then 
          jo_lagged=1 
      else 
          jo_lagged=2
      end if 
      
      

     if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
        

    
do jzt=1, n_zt
    if ((jzt.ne.1).AND.(jo.eq.1)) exit
     
!************************************************************************************************************
!     CALCULATE IF YOU ARE ELIGIBLE FOR TRANSFER AND MAX SAVINGS, AND YOUR INCOME GROUP
!************************************************************************************************************
!initialize
  actual_earn =  0.0d0
  income_bt=0.0d0
  tax_SS=0.0d0
  tax_med=0.0d0
  income_pay(:)=0.0d0
  Total_tax_pay(:)=0.0d0
  base_2_pay(:)=0.0d0
  base_1_pay(:)=0.0d0
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay(:)=0.0d0
  V_c_all_trpay(:)=0.0d0
  m=0.0d0
  n=0.0d0
  SAVE_EXPECTED_B_PRIME_trpay(:,:)=0.0d0 ! this changes with every new state considered. but depends on jsav.
SAVE_EXP_L_prime_INT(:,:,:,:,:,:)=0.0d0
  TR_indic=0
  

! earnings and income     
  if (jo.eq.1)  then 
      earnings_interpol=0.0d0
  else 
      call InterpolateScalar ( HK(0 : n_hc), base_earn(0:n_hc,jh,je,jfix,jzt,jo) , xfn(jhc,jage)  , earnings_interpol ) 
  end if 

  
  actual_earn = max(0.0d0 , ( earnings_interpol * (hrs_offer(jo) -phi(je,jh,jtype2) )) ) ! earnings net of sick days
  income_bt= ((interest-1.0d0))*as(jas,jtype1,jage) +  actual_earn  ! income before tax is equal to interest income plus earnings net of sick days
  
  if ((jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23)) THEN
      
        if (jo.ne.1) then ! husband works - deduct from his income
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  ))
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife. 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)))
end if 
      
      else 
  if ((jo.eq.3).OR.(jo.eq.5)) then ! husband holds insurance -- so deduct from his income only
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife. note this can be zero
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)))
end if 
  
end if 
      
  
  ! determine income group
if (income_bt.le.INC_TS(1)) then 
jinc=1
else if (income_bt.le.INC_TS(2)) then 
jinc=2
else if (income_bt.le.INC_TS(3)) then 
jinc=3
else if (income_bt.le.INC_TS(4)) then 
jinc=4
else 
jinc=5
end if

!income and taxes assuming you pay. we assume that you always pay this mimimum costs if no shocks
  income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo) + oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))  +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(1) =  prem_ephi(jo,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo) +  oop_spouse(jage,je,jtype4,jo) +Total_tax_pay(1)  ! total resources subtracted 
  base_1_pay(1) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        ! cash on hand (if you pay the OOP)
  
! solve for the EV associated with having 0 assets, and getting treated, and also 0 assets but not getting treated.
  EXPECTED_B_PRIME_C_FLOOR=0.0d0
  EXPECTED_B_PRIME_C_FLOOR_NT =0.0d0 
  EXP_L_prime_INT=0.0d0
  check=0.0d0

        do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
               call TYPE3_DET(jh_p,jr_p,jtype3_p) 
            !INTERPOLATE FOR HUMAN CAPITAL
            if (HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1) .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(1,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(1,jtype3_p, (HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(1,jtype3_p, HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1), jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1))+1, jage+1)  -   xfn((HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)), jage+1))
		    inter_hc = EXP_L_prime_simple(1,jtype3_p, HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1), jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
             

            EXPECTED_B_PRIME_C_FLOOR= EXPECTED_B_PRIME_C_FLOOR      + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) * (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * (cost_death + Bequest_val(0.0d0) )  ) 
            EXPECTED_B_PRIME_C_FLOOR_NT=EXPECTED_B_PRIME_C_FLOOR_NT + prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) * (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * (cost_death + Bequest_val(0.0d0) )  ) 
        end do ; end do ; end do; end do ; end do    

        
                 
!Calculate value if on cons floor - here you get treated  
if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then
 
    CON=cons_floor(jage, jmar,je,jh)
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON,UTI)   ! the utility is that associated with paying and being treated. 
    B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = UTI+ beta_hat(je)*EXPECTED_B_PRIME_C_FLOOR
           

else   

    MAX_SAV_pay(1) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  -cons_floor_at(jage, jmar,je,jh) 

          jsav_upper_pay(1)=0
          do while ( MAX_SAV_pay(1) .ge.  sav(jo,jsav_upper_pay(1)+1))
           jsav_upper_pay(1)=jsav_upper_pay(1) +1
           if (jsav_upper_pay(1).EQ.n_sav) exit
          end do 
          
    
 ! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) 
     V_cf_trpay(1)=0.0d0
     ASP= MAX_SAV_pay(1) +as(jas,jtype1,jage)
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)      
    !Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0   
 EXP_L_prime_INT=0.0d0
     

      if (ASP .lt. as(AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1)),jtype1,jage+1)) then 
          pause 
      else  
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))
          
          do while (ASP .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      end if 
      
     do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)  
          

        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)
    ! m was set outside the loops
      

        if (as(n_as,jtype1,jage+1) <= ASP) then
            !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn(n, jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 78 
 
        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 78             
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 78                
  
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) + &
                mx1*(1-mx2)*yt(2) + &
                (1-mx1)*mx2*yt(3) + &
                mx1*mx2*yt(4) 

78    EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo)  *(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death )
    end do; end do  ; end do;   end do; end do
      V_cf_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(1)

! Calculate value of consuming all and saving nothing
     V_c_all_trpay(1) =0.0d0
     call func_uti(jage,1,je,jh,jo,jtype4,jtype2, base_1_pay(1)/(1.0d0+taxc), UTI)  
     V_c_all_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR ! this was solved earlier as the EV of saving nothing and getting treated    
    


! BEGIN SAVINGS LOOP 
 maxB_trpay(1)=V_cf_trpay(1)
do jsav= jsav_upper_pay(1) , 1, -1 
                 
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
        if (ASP.le.0.d0)  exit ! we already did the case where consume all and save 0
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
    CON_pay(1)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav(jo,jsav) )/(1.0d0+taxc)
    

    
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON_pay(1),UTI)

!Find expected value
 EXPECTED_B_PRIME_trpay(:)=0.d0   
 EXP_L_prime_INT=0.0d0  
 VF=0.0d0
      do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)                    
          
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        m=AS_LB(jas,jtype1,jage+1,jo,jsav)
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1) 
        
        

        if (as(n_as,jtype1,jage+1) <= ASP) then
            !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn(n, jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 79 

        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 79             
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 79                
  
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) + &
                mx1*(1-mx2)*yt(2) + &
                (1-mx1)*mx2*yt(3) + &
                mx1*mx2*yt(4) 

        
79        SAVE_EXP_L_prime_INT(jsav,jh_p,jr_p, jo_p,jmar_p,jzhc_p) =  EXP_L_prime_INT
         

            EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)    
      end do ; end do ; end do   ; end do ; end do 

SAVE_EXPECTED_B_PRIME_trpay(jsav,1)=EXPECTED_B_PRIME_trpay(1) ! we save this EV of getting treated for every jas. we will use it later when solving for value of jtrpay=2 (treat not pay)
     
VF=UTI + beta_hat(je)*EXPECTED_B_PRIME_trpay(1)


     if (VF.lt.maxB_trpay(1)) exit
        
        if ((VF).ge.(maxB_trpay(1))) then                              
              maxB_trpay(1)=VF  
        end if  
                                
end do  !jsav 
    
! now find the maximum considering the extremes 
 
        B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(1),  V_c_all_trpay(1)) 

   
end if   ! for being on cons floor     
    
 !if no shocks, we only need to solve the value if you "pay and get treated"
if (jtype2.eq.1) then      
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 
 ! *******************************************************************************************************************************************************   
else 
!IF YOU HAVE HEALTH SHOCK, YOU HAVE TO COMPARE not TREAT AND not PAY OPTIONS
!income and taxes if not pay
    
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))  +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(2) =  prem_ephi(jo,jtype4) +  oop_spouse(jage,je,jtype4,jo) + Total_tax_pay(2)  ! total resources subtracted 
  base_1_pay(2) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        ! cash on hand (if you do not pay the OOP)


      
!Calculate value if not on cons floor    
 if ( base_1_pay(2) .gt. (cons_floor_at(jage, jmar,je,jh))) then 
 
!Find max savings if you don't pay
   MAX_SAV_pay(2) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(2)  -cons_floor_at(jage, jmar,je,jh) ! if not pay OOP
        
          !find lower bound 
          jsav_upper_pay(2)=0
          do while ( MAX_SAV_pay(2) .ge.  sav(jo,jsav_upper_pay(2)+1))
           jsav_upper_pay(2)=jsav_upper_pay(2) +1
           if (jsav_upper_pay(2).EQ.n_sav) exit
          end do 
          
      
! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) - here, you don't pay oop and get treated; and also solve when  you don't get treated.
     V_cf_trpay(2)=0.0d0
     V_cf_trpay(3)=0.0d0
     ASP= MAX_SAV_pay(2) +as(jas,jtype1,jage)  
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
     !Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0   
 EXP_L_prime_INT=0.0d0
   
     
      if (ASP .lt. as(AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2)),jtype1,jage+1)) then 
          pause 
      else  
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2))
          
          do while (ASP .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      end if 
      
 
     do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)  
                  
          
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)
        
        if (ASP ==0.0d0 .and. HK_INTERPOL(jzhc_p,jo) == 0.0d0) then
           EXP_L_prime_INT= EXP_L_prime_simple(1,jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)
           go to 87
        end if
        if (ASP ==0.0d0) then	
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(1,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(1,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(1,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(1,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if      
            go to 87             
        end if
        if (as(n_as,jtype1,jage+1) <= ASP) then
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 87 

        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 87             
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 87                
 
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) + &
                mx1*(1-mx2)*yt(2) + &
                (1-mx1)*mx2*yt(3) + &
                mx1*mx2*yt(4) 
 
            
87    EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
    EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) + prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo)*(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)

     end do; end do  ; end do;   end do; end do
         call func_uti(jage,2,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_cf_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(2)  
         call func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_cf_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(3)
  
! Calculate value of consuming all and saving nothing - don't pay, get treated; and also for don't get treated
     V_c_all_trpay(2) =0.0d0
     call func_uti(jage,2,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI)  
     V_c_all_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR ! this was solved earlier as the EV of saving nothing and getting treated  
  
  !if not get treated
       V_c_all_trpay(3) =0.0d0
     call func_uti(jage,3,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI) 
  V_c_all_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 
  
   
! BEGIN SAVINGS LOOP   
  !we initiate the max to be the value associated with saving all you can save, ie, V_cf_trpay()
  maxB_trpay(2)=V_cf_trpay(2)
  maxB_trpay(3)=V_cf_trpay(3)
do jsav= jsav_upper_pay(2) , 1, -1 
                    
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
     if (ASP.le.0.d0)  exit  
    CON_pay(2)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(2) -sav(jo,jsav) )/(1.0d0+taxc)
     


                    call  func_uti(jage,2,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(2) = UTI
                    call  func_uti(jage,3,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(3) = UTI   
                
 VF_tr(2:3)=0.0d0                                        
!CALCULATE EXPECTED VALUE 
if ((SAVE_EXPECTED_B_PRIME_trpay(jsav,1)).ne.(0.0d0)) then ! we already solved for the EV of treating. and we can use the saved interpolated values to find value of not treating too.
    VF_tr(2)=UTI_trpay(2) + beta_hat(je)*SAVE_EXPECTED_B_PRIME_trpay(jsav,1)
  
! solve EV for not pay not treat 
    EXPECTED_B_PRIME_trpay(3)=0.0d0
       do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)    
        EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) +  prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* SAVE_EXP_L_prime_INT(jsav,jh_p,jr_p, jo_p,jmar_p,jzhc_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
        end do; end do  ; end do;   end do; end do    
            
    VF_tr(3)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
                
                
                
else 
 ! for jsav we have not considered before, solve EV               
            EXP_L_prime_INT=0.0d0
            EXPECTED_B_PRIME_trpay(:)=0.0d0   
             m=AS_LB(jas,jtype1,jage+1,jo,jsav)
             
     do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)                    
          
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)

        if (as(n_as,jtype1,jage+1) <= ASP) then
            !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 88 
        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 88            
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 88                
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) +  mx1*(1-mx2)*yt(2) + (1-mx1)*mx2*yt(3) +  mx1*mx2*yt(4) 

            
88            prob_d=0.0d0
            prob_d=(1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death
            
            EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  prob_d)
            EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) +  prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  prob_d)

     end do ; end do ; end do   ; end do ; end do 
     

           VF_tr(2)=UTI_trpay(2) + beta_hat(je)*EXPECTED_B_PRIME_trpay(2)
           VF_tr(3)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
           
end if 



! if both values are going down, stop and exit
     if (((VF_tr(2)).lt.(maxB_trpay(2))).AND.(VF_tr(3).LT.maxB_trpay(3)))  exit  
         
              if (VF_tr(2).ge.maxB_trpay(2)) then                              
              maxB_trpay(2)=VF_tr(2)   
              end if   
              
              if (VF_tr(3).ge.maxB_trpay(3)) then                              
              maxB_trpay(3)=VF_tr(3)  ! will save the new max and continue iterating over savings. 
              end if 
                      
           
         
end do   !for jsav    

   
    
    ! now find the maximum considering the extremes 
  
        B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(2),  V_c_all_trpay(2)) 
        B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(3),  V_c_all_trpay(3)) 

    
               
 
 else 
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 !  solve value of being on cons floor and not treating
          call func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 

          
 end if   ! for not being on cons floor
 
 if ((I_treat_Mcaid.EQ.1).AND.( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) ) THEN
      B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = max( B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) , B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)) 
 END IF 
 
end if  

if (jtype2.eq.1) then
 B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) ) 
else 
    if ((Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage)).NE.0.0d0) then 
B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = (Shock_type(2,jo,jh,jage)/(Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage)))* max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) ) + (Shock_type(1,jo,jh,jage)/(Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage))) *  B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)
end if
end if 

  

if (jExperimentNumber.eq.13) THEN ! force everyone to treat and pay
    B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) =  B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)
    else 
B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage), B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) ,B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) )
    end if 
    
end do ! jzt
end do ! jo
B_1(1,jMEcat,jtype4,2:n_zt,jas,jtype3,jtype2,jhc,jage) = B_1(1,jMEcat,jtype4,1,jas,jtype3,jtype2,jhc,jage) !for those who don't get an offer, zt does not matter, so set all equal to when zt =1
B_2(1,jMEcat,jtype4,2:n_zt,jas,jtype3,jtype2,jhc,jage) = B_2(1,jMEcat,jtype4,1,jas,jtype3,jtype2,jhc,jage) !for those who don't get an offer, zt does not matter, so set all equal to when zt =1

end do ! jMEcat
end do ! jmar
end do ! jas
end do ! jr
end do; end do; end do ! for  du, dp, s


do jas =1, n_as ; do jr=1, n_r ;  do jzt=1, n_zt ;  do jtype4=1, n_type4
 call TYPE3_DET(jh,jr,jtype3)  ! Gives Type3 given jh and jr

! Find "EXPECTED_B" over health shocks, ME_cat
        EXPECTED_B(:)=0.0d0 
    do jo=1, n_o
            do jtype2=1, n_type2 
            do jMEcat=1, n_MEcat
                       
                prob_s=0.0d0
                prob_s= shocks_risk(je,jtype2,jh,jr, jage)*Prob_MEcat(jMEcat)
   EXPECTED_B(jo)= EXPECTED_B(jo)+  prob_s * (   (Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage)) * B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) + Shock_type(3,jo,jh,jage) * B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  )
            end do 
            end do 
            
            B(jo,jtype4,jzt,jas,jtype3,jtype1,jhc,jage) = EXPECTED_B(jo) 

    
    end do  
              
           
 ! decide if accept offer and save EXP_L_prime_fix
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,:,jtype4)=EXPECTED_B(1) ! initialize to the non-employed value
    
    If   (EXPECTED_B(5).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,5,jtype4)=EXPECTED_B(5)
    end if
  
    If   (EXPECTED_B(4).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,4,jtype4)=EXPECTED_B(4)
    end if  
    
    If   (EXPECTED_B(3).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,3,jtype4)=EXPECTED_B(3)
    end if
  
    If   (EXPECTED_B(2).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,2,jtype4)=EXPECTED_B(2)
    end if 

    
end do ; end do; end do ; end do  
end do 
end do hkloop


!Get "EXP_L_prime_simple" which does not depend on jos (wife working) and jzt. We take probabilities over these. 
do jo=1,n_o  ;  do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid   
! get back jh here because it's needed for probability of wife working
        if ((jtype3.eq.1).OR.(jtype3.eq.2)) then 
        jh=1
        else if ((jtype3.eq.3).OR.(jtype3.eq.4)) then
    jh=2
    else if ((jtype3.eq.5).OR.(jtype3.eq.6)) then
        jh=3
    end if 
    
    ! NOT EMPLOYED LAST YEAR
   EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,:,1) =0.0d0
    do jzt=1,n_zt ! for singles and married
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,1) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,1)  + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,1) * p_zt(je,jzt,1)
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,1) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,1)  +   p_zt(je,jzt,1)* (EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,2)*prob_os(jage,je,jh,1,jfix)   + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,3)*prob_os(jage,je,jh,2,jfix) ) 
    end do
   
    !EMPLOYED LAST YEAR
       EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,:,2) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,2) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,2)  + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,1) * p_zt(je,jzt,2)
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,2) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,2)  +   p_zt(je,jzt,2)* (EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,2)*prob_os(jage,je,jh,1,jfix)   + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,3)*prob_os(jage,je,jh,2,jfix) ) 
    end do
    
end do; end do; end do; end do 

end do ageloop 



do jo=1,n_o  ; do jmar =1,n_mar ; do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid  ; do jage=1,n_ret-2
EXP_L_prime(jas,jtype1,jtype3,jhc,jage,jo,jmar,1) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,jmar,1)
EXP_L_prime(jas,jtype1,jtype3,jhc,jage,jo,jmar,2) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,jmar,2)
end do; end do; end do; end do; end do; end do


end do ! jht
end do ! jfix



!***********************************************************************
! SIMULATE ECONOMY,
 Call SIMULATION(EXP_L_prime, B,  RET_PRIME ,  jExperimentNumber )
!**************************************************************

end do ! for experiments

End Program VALUEFUNCTION